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Abstract:  This  report  describes  a  particle  tracking  computer  program 
named  PT123.  The  development  of  PT123  was  supported  in  part  by  the 
Civil  Works  Basic  Research  project  entitled  “Efficient  Resolution  of 
Complex  Transport  Phenomena  Using  Eulerian-Lagrangian  Techniques” 
and  in  part  by  the  System- Wide  Water  Resources  Program  (SWWRP). 
Given  velocities,  PT123  can  track  massless  particles  in  1-,  2-,  and  3-D 
unstructured  or  converted  structured  meshes.  The  elements  used  to  con¬ 
struct  PT123  meshes  are  line  elements  in  l-D,  triangular  and/or  quadri¬ 
lateral  elements  in  2-D,  and  tetrahedral,  triangular  prism,  and/or  hexa- 
hedral  elements  in  3-D.  One  adaptive  (embedded  4th-  and  5th-order)  and 
three  non-adaptive  (1st-,  2nd-,  and  4th-order)  Runge-Kutta  (RK)  methods 
are  included  in  PT123  to  solve  the  ordinary  differential  equations  describ¬ 
ing  the  motion  of  particles.  The  adaptive  RK  method  allows  the  user  to 
control  tracking  accuracy  with  specified  error  tolerances.  The  non- 
adaptive  RK  methods  provide  the  user  options  to  balance  computational 
efficiency  and  accuracy  by  using  lower  order  schemes  for  smooth  velocity 
fields  and  higher  order  schemes  for  complex  velocity  fields.  Both  element- 
by-element  (EBE)  and  non-element-by-element  (NEBE)  tracking 
approaches  are  incorporated  into  PT123.  Both  node-  and  element-based 
velocity  can  be  used  for  particle  tracking.  PT123  can  execute  forward  and 
backward  tracking  and  output  tracking  history  at  a  specified  frequency.  It 
tracks  particles  along  the  closed  boundary  and  stops  tracking  when  a 
particle  encounters  the  open  boundary  through  which  particles  enter  or 
exit  the  computational  domain.  The  start  and  end  times  of  tracking  are 
flexible  as  long  as  their  corresponding  velocities  can  be  computed  via 
temporal  interpolation  using  the  given  velocities.  This  report  is  the  first 
report  of  the  series  describing  the  development  and  application  of  PT123. 
It  details  the  governing  equation  and  numerical  approaching  associated 
with  PT123  Version  1.0.  Six  test  examples  in  multiple  dimensions  are  used 
for  verification  and  demonstration.  The  structure  and  the  input  guide  of 
the  computer  program  are  given  in  the  appendices. 
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1  Introduction 

The  particle  tracking  (PT)  technique  has  a  wide  range  of  application  in 
environmental  sciences  and  engineering.  This  technique  typically  uses  the 
output  from  hydrodynamic  and/or  advection-diffusion  models  to  predict 
particle  movements  in  a  Lagrangian  manner.  Given  the  velocity  field,  PT 
can  provide  a  quick  estimate  of  how  a  chemical  migrates  in  complex  sur¬ 
face  water  and  groundwater  systems.  It  can  be  used  to  understand, 
visualize,  and  analyze  flow  fields  (Pokrajac  and  Lazic  2002).  It  can  be  used 
to  study  sediment  transport  (MacDonald  et  al.  2006),  oil  spill  (Liu  et  al. 
2011),  and  natural  or  man-induced  retardation  mechanisms  that  maybe 
used  for  the  remediation  or  prevention  of  environmental  pollution.  It  can 
be  used  to  understand  and  predict  fish  behavior  (Goodwin  et  al.  2006)  for 
ecosystem  restoration  and  preservation.  It  can  also  be  applied  in  the 
Eulerian-Lagrangian  (EL)  approximation  to  solve  transport  equations 
numerically,  which  is  a  crucial  modeling  practice  to  help  deal  with  envi¬ 
ronmental  issues  concerning  water  quality.  The  quality  of  particle  tracking 
dictates  much  of  the  accuracy  of  the  whole  EL  approximation  (Russell  and 
Celia  2002)  as  well  as  efficiency  on  serial  and  parallel  platforms  (Cheng 
and  Plassman  2004).  Efficient  numerical  methods  for  transport  are  neces¬ 
sary  for  large-scale  modeling  in  achieving  the  Corps’  mission.  Therefore, 
having  accurate  and  efficient  PT  can  help  the  Corps’  engineers  and  scien¬ 
tists  to  carry  out  reimbursable  and  research  and  development  (R&D) 
projects.  These  methods  have  been  implemented  in  a  computer  model 
called  PT123. 

1.1  Purposes  of  PT123  research  study 

The  purposes  of  PT123  research  are  two-fold.  One  is  to  construct  accurate 
and  efficient  PT  computer  routines  for  solving  multi-dimensional  trans¬ 
port  problem  using  the  Eulerian-Lagrangian  localized  adjoint  methods 
(ELLAM)  numerical  method  (Russell  and  Celia  2002),  as  proposed  in  the 
Civil  Works  Basic  Research  project  entitled  “Efficient  Resolution  of 
Complex  Transport  Phenomena  Using  Eulerian-Lagrangian  Techniques.” 
The  other  is  to  develop  a  library-type  computer  program  which  can  be 
incorporated  easily  into  or  linked  to  ERDC’s  existing  flow  or  transport 
models  (e.g.,  adaptive  hydraulics  model  (ADH),  particle  tracking  model 
(PTM))  to  enhance  computational  accuracy  and  efficiency  in  various 
applications. 
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1.2  Modeling  approach 

The  velocity  field  dictates  the  PT  result.  Inaccuracies  in  velocity  values 
introduce  an  error  into  PT.  Spatial  interpolation  of  velocities  introduces 
another  error  into  PT  because  the  exact  velocity  field  differs  from  the 
interpolated  field  even  if  the  nodal  velocities  are  exact  (Pokrajac  and  Lazic 
2002). 

Analytical  PT  solutions  are  limited  to  the  cases  with  simple  geometry  and 
velocity  fields.  Semi-analytical  PT  is  used  in  the  Pollock’s  method,  where 
linear  interpolation  of  velocity  enables  the  analytical  calculation  of  path 
lines  and  travel  times  over  an  element  (Pollock  1988).  The  standard 
numerical  methods  used  are  the  ist-order  Euler,  the  2nd-order  Euler 
predictor-corrector,  or  higher-order  Runge-Kutta  (RK)  methods.  For 
convenience,  RKi  and  RK2  are  used  to  represent  the  1st-  and  2nd-order 
methods  (Press  et  al.  1992).  The  previous  studies  showed  that  higher- 
order  RK  methods,  e.g.,  4st-order  RK  or  RK4,  are  superior  to  the  lower- 
order  RK  methods  regarding  accuracy,  and  adaptive  spatial  or  temporal 
steps  improve  significantly  the  efficiency  of  PT  algorithms  in  velocity  fields 
containing  wide  spectrums  of  velocity  magnitude  and  element  size  (Cash 
1989;  Press  et  al.  1992;  Bensabat  et  al.  1998;  Oliveira  and  Baptista  1998); 
Cheng  and  Plassman  2004. 

By  assuming  that  the  given  velocity  field  is  accurate  and  the  velocity  inter¬ 
polation  error  is  negligible,  the  PT123  implementation  focuses  primarily 
on  the  techniques  for  solving  the  ordinary  differential  equations  (ODEs) 
that  describe  the  motion  of  particles  along  their  path  lines.  The  PT123 
computer  program  presented  in  this  document  is  the  result  of  the  initial 
effort  of  PT123  basic  R&D,  resulting  in  a  stand-alone  computer  code.  The 
next  planned  research  tasks  for  PT123  include 

1.  parallelization; 

2.  GUI  development; 

3.  library  format; 

4.  incorporation  of  mechanisms/processes  that  modify  tracking 
velocities; 

5.  conversion  of  finite  difference  or  finite  volume  model  data  into  PT123 
format. 
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Refined  and  tested  computational  algorithms  of  PT123  will  be  ported  into 
ERDC’s  in-house  models,  e.g.,  ADH  (ADH  2010)  and  PTM  (MacDonald 
et  al.  2006).  A  summary  of  PT123  computational  strategy  and  feature  and 
model  input/output  (I/O)  follows  next. 

1.3  Computational  strategy  and  features 

PT123  employs  the  following  techniques  in  its  computation: 

1.  Perform  PT  on  either  an  element-by-element  (EBE)  or  a  non-element- 
by-element  (NEBE)  basis. 

2.  Use  both  absolute  tolerance  (ATOL)  and  relative  tolerance  (RTOL)  to 
control  accuracy  in  time  integration  when  adaptive  RK  is  used. 

3.  Use  interpolation  to  estimate  the  derivative  term,  i.e.,  velocity,  during 
the  RK  process.  The  interpolation  is  linear  in  time  and  consistent  with 
element  shape  function  in  space. 

With  the  strategy  implemented,  the  current  version  (1.0)  of  PT123  includes 
the  following  computational  features: 

1.  PT  in  multiple  dimensions:  1-,  2-,  or  3-D. 

2.  Flexible  time  and  length  units:  any  time  and  length  units  are  valid  if 
consistent. 

3.  Different  RK  schemes: 

a.  One  adaptive:  embedded  4th-  and  5th-order  (RK45). 

b.  Three  non-adaptive:  ist-order  (RKi),  2nd-order  (RK2),  and 
4th-order  (RK4)  schemes  can  be  used  at  the  user’s  choice. 

4.  Forward  or  backward  PT :  the  user  specifies  in  the  input  data  whether 
forward  or  backward  PT  is  to  be  performed. 

5.  Multiple  particles:  there  is  no  constraint  on  the  number  of  tracked  par¬ 
ticles;  the  user  specifies  the  number  and  the  locations  of  the  tracked 
particles  in  the  input  data. 

6.  Steady  or  transient  velocity  fields:  the  user  specifies  whether  a  steady 
or  transient  velocity  is  to  be  applied  for  PT. 

7.  Node-  or  element -based  velocity:  when  node-based  velocity  is  used, 
only  one  velocity  is  assigned  to  a  global  node  at  a  time;  when  element- 
based  velocity  is  used,  the  velocity  at  a  node  may  change  with  the 
element  involving  that  node  due  to  heterogeneity  or  other  reasons. 

8.  Velocity  conversion  factor:  each  node  or  each  element  is  assigned  a 
velocity  conversion  factor  to  allow  the  conversion  from  the  given  flow 
velocity  to  the  tracking  velocity,  e.g.,  from  the  given  Darcy  velocity  to 
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the  pore  velocity  for  tracking.  This  feature  offers  the  flexibility  of  using 
various  tracking  velocities  for  particles  of  different  kinds,  e.g.,  sedi¬ 
ments  of  different  sizes,  in  the  same  flow  field. 

9.  Courant  number  to  control  PT  time  step  size:  a  user-specified  Courant 
number  value  can  be  used  to  compute  the  PT  time  step  size  using  the 
tracking  velocities  and  the  characteristic  length  associated  with  the 
element  that  contains  the  particle  being  tracked. 

10.  Flexible  start  time  and  time  duration  for  tracking:  the  time  parameters 
can  be  assigned  any  values  in  the  range  that  the  given  velocity  field 
covers. 

11.  Various  types  of  element  shapes:  PT123  can  compute  PT  within 
unstructured  meshes  composed  of  line  elements  in  l-D,  triangular  and 
quadrilateral  elements  in  2-D,  and  tetrahedral,  triangular  prism,  and 
hexahedral  elements  in  3-D  domains. 

1.4  Input  and  output 

PT123  does  not  require  any  specification  of  time  and  length  units  in  the 
input  file.  Any  combination  of  time  and  length  units  can  be  utilized  in  PT 
computation  as  long  as  consistency  is  maintained  throughout  the  input 
data.  The  output  assumes  the  same  time  and  length  units  implied  in  the 
input  data. 

The  input  data  that  PT123  requires  for  PT  computation  includes 

1.  element  indices  and  nodal  coordinates  (i.e.,  geometry  of  the 
computational  domain); 

2.  velocities  (e.g.,  flow  fields  from  hydrodynamic  models); 

3.  velocity  conversion  factors; 

4.  open/closed  boundary  information; 

5.  specifics  for  PT  computation  (e.g.,  particle  data,  computation 
parameters,  etc.). 

PT123  uses  several  input  files  to  accommodate  its  input  data.  The  details 
of  the  PT123  input  files  are  given  in  Appendix  B. 

PT123  outputs  the  trajectory  of  each  tracked  particle  from  the  start  time 
through  the  time  duration,  i.e.,  time  versus  location  for  each  particle,  at  a 
desired  frequency.  For  example,  if  a  particle  is  tracked  for  a  time  duration 
of  30,000  sec  and  the  user  wants  to  trace  the  locations  of  the  tracked 
particle  every  100  sec,  then  the  trajectory  will  include  301  points,  where 


ERDC  TR-11-10 


5 


301  is  equal  to  (30,ooo/ioo)+i.  PT123  stores  the  PT  output  in  ASCII  and 
BINARY  format  for  inspection  and  post-processing,  respectively. 

The  remainder  of  the  report  provides  details  of  the  information  summa¬ 
rized  so  far.  The  mathematical  statements  and  numerical  solutions  are 
stated  in  Chapter  2.  Six  test  examples  are  given  in  Chapter  3  for  both 
verification  and  demonstration  purposes.  Final  remarks  on  the  develop¬ 
ment  of  PT123  and  an  outline  of  tasks  for  future  advancements  are  given 
in  Chapter  4.  The  program  structure  and  subroutine  description  are  pro¬ 
vide  in  Appendix  A.  The  input  guide  is  given  in  Appendix  B.  The  output 
files  are  described  in  Appendix  C. 
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2  Governing  Equations  and  Numerical 
Solutions 

2.1  Governing  equation 

In  PT123,  the  following  ODE  in  vector  form  is  solved  for  defining  the 
particle  path. 

^  =  V(x,t)  (1) 

at 

where: 

d  =  Courant  number 
x  =  location  of  a  tracked  particle  [L] 
t  =  time  [t] 

V  =  tracking  velocity  [L/t]. 

Given  the  initial  location  of  a  particle,  i.e.,  x(fo),  one  can  compute  the 
particle  path  through  time  integration  over  the  specified  velocities,  as 
shown  in  Equation  2. 

t 

x(t)  =  x(t0)  +  J  V(x,t')dt'  (2) 

*0 

where: 

to  =  start  time  for  PT  [t] 

t'  =  a  dummy  variable  used  for  time  integration. 

2.2  Time  integration 

PT123  includes  an  adaptive  time  integration  algorithm,  where  the  differ¬ 
ence  from  embedded  4th-  and  5th-order  RK  results  is  employed  for  error 
estimation.  The  estimated  error  is  compared  to  the  prescribed  error 
tolerances  to  adjust  the  time  step  size  for  time  integration  in  the  PT  com¬ 
putation.  PT123  also  provides  options  of  using  1st-,  2nd-,  or  4th-order  RK 
for  computation,  where  the  user  provides  a  specified  time  step  size,  i.e., 
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the  DT_INITo  parameter  in  the  PT  Specifics  file,  in  Appendix  B.  The  com¬ 
puted  particle  trajectory  is  composed  of  many  tracking  segments,  and  each 
segment  is  associated  with  a  successful  RK  computation. 

2.2.1  Adaptive  RK  schemes 

With  RK45,  we  first  compute  k,  (i  =  1  to  6)  as  defined  as  in  Equation  3, 
where  k,  is  a  vector  of  a  size  equal  to  the  dimension  of  the  PT  spatial 
domain,  i.e.,  1  for  l-D,  2  for  2-D,  and  3  for  3-D.  Table  1  lists  the  values  of 
coefficients  a,-,  6,7,  a,  and  caused  in  RK45  as  shown  below  (Press  et  al. 

1992). 

Table  1.  Cash-Karp  coefficients  for  the  embedded  4th-  and  5th-order  RK  (from  Press  et  al.  1992). 


/ 

3/ 

bj 

Ci 

C* 

1 

0 

0 

0 

0 

0 

0 

37/387 

2825/27648 

2 

1/5 

1/5 

0 

0 

0 

0 

0 

0 

3 

3/10 

3/40 

9/40 

0 

0 

0 

250/621 

18575/48384 

4 

3/5 

3/10 

-9/10 

6/5 

0 

0 

125/594 

13525/55296 

5 

1 

-11/54 

5/2 

-70/27 

35/27 

0 

0 

277/14336 

6 

7/8 

1631/55296 

175/512 

575/13824 

44275/110592 

253/4096 

512/1771 

1/4 

j  = 

1 

2 

3 

4 

5 

k1=Af-V(xn,fJ 

k2  =  Af  •  V(xn  +  b21k1,tn  +  a2  Af) 

k3  =  At-V(xn  +  b31k1  +  b32k2,tn  +  a3At) 

k4  =  Af  •  V(xn  +  h41k1  +  642k2  +  b43k3 ,  tn  +  a4  Af ) 

k5  =  Af  ■  V(xn  +  h51k1  +  h52k2  +  h53k3  +  b54k4 ,  tn  +  a5  Af ) 

k6  =  Af  V(xn  +h61  k4  +h62k2  +h63k3  +h64k4  +b65k5,fn  +a6Af) 


where: 

Af  =  step  size  for  time  integration 
xn  =  start  location. 

Therefore,  the  embedded  5th-order  RK  yields 

Xn+1  =  Xn  +  Clkl  +  C2k2  +  C3k3  +  C4k4  +  C5k5  +  C6k6  +  0(Af  ®)  (4) 
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While  the  embedded  4th-order  results  in 

XL+1  =  Xn  +  <hkl  +  C2k2  +  C3k3  +  <bk4  +  C5k5  +  C6k6  +  0(Af5)  (5) 


where: 

xn  i  =  estimated  end  location  using  the  embedded  5th-order  RK 
xn  fl  =  estimated  end  location  using  the  embedded  4th-order  RK. 

2.2.2  Error  estimate 

Using  RK45,  the  integration  error  can  be  estimated  using  Equation  6. 


A  = 


Xn+1  Xn+1 


(6) 


where: 


A  =  estimated  error  of  time  integration. 


2.2.3  Adaption  of  time  step  size 

Equation  7  is  used  to  compare  the  estimated  error  with  prescribed  error 
tolerances  for  the  adaption  of  time  step  size. 


<5,= 


ATOL  +  RTOL  ■  max 


Xn+lJ  XnJ 


Xn+lJ-XnJ 


(7) 


where: 

8j  =  the  ratio  of  the  estimated  error  to  the  prescribed  error 
tolerance  in  the  /-th  spatial  direction 
ATOL  =  prescribed  absolute  error  tolerance  [L] 

RTOL  =  prescribed  relative  error  tolerance  [dimensionless] 
xn+ij  =  j-th  component  of  xn+1 . 


Here  ATOL  represents  the  allowed  estimated  error  of  time  integration  for 
each  tracking  segment  in  an  absolute  sense.  On  the  other  hand,  RTOL  is 
the  allowed  error  portion  when  compared  to  the  length  of  the  tracking 
segment  being  calculated.  The  combination  of  the  two,  as  described  in  the 
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denominator  on  the  right-hand  side  of  Equation  7,  defines  the  allowed 
error  tolerance  for  each  tracking  segment.  The  user  chooses  ATOL  and 
RTOL  based  on  the  requirement  of  accuracy  for  his/her  application. 

We  now  define  a  ratio  as 


R  = 


1 

mays,.) 


(8) 


and  use  R  in  determining  an  appropriate  time  step  size  for  PT  compu¬ 
tation.  Two  possibilities  exist: 

When  R  <  1: 

When  R  is  smaller  than  1,  the  estimated  error  exceeds  the  allowed  error 
tolerance,  i.e.,  Equation  7.  In  this  case,  the  time  step  size  is  reduced  using 
the  following  equation: 


At*  =A  t-SFR025 


(9) 


where: 

At'  =  adapted  PT  time  step  size 
SF  =  safety  factor  used  in  adaption. 

When  R  >  1: 

When  R  is  greater  than  or  equal  to  1,  the  estimated  error  is  smaller  than  or 
equal  to  the  allowed  error  tolerance.  In  this  case,  the  time  step  size  for  the 
current  PT  computation  is  small  enough  to  meet  the  required  accuracy, 
and  the  time  step  size  can  be  increased  for  the  successive  particle  tracking 
computation.  This  increased  time  step  size  is  computed  using  the  follow¬ 
ing  equation: 


At*  =AtSFR02 


(10) 


2.3  Interpolation  of  velocity 

The  given  velocity  field  can  vary  in  both  time  and  space.  While  the  ana¬ 
lytical  velocity  is  not  available  in  complex  real-world  systems,  interpola- 
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tion  becomes  essential  to  estimate  velocity  at  various  times  and  locations 
in  the  PT  computation.  PT123  conducts  a  linear  temporal  interpolation 
and  the  following  spatial  interpolation  schemes,  depending  on  the  shape  of 
the  active  element,  for  its  velocity  computation: 

1.  Linear  for  l-D  line,  2-D  triangular,  and  3-D  tetrahedral  elements 

2.  Bi-linear  for  2-D  quadrilateral  elements 

3.  tri-linear  for  3-D  hexahedral  elements 

4.  Combined  linear /bi-linear  for  3-D  triangular  prism  elements. 

2.4  Element-by-Element  (EBE)  tracking 

PT123  can  conduct  PT  on  an  EBE  basis  (Cheng  et  al.  1996),  where  each 
tracking  segment  computed  using  the  designated  RK  scheme  is  within  an 
element.  Figure  1  presents  this  EBE  tracking  concept. 


Tracking  ^ 
Time  Consumed 
v^CompletelyJ^ 


Particle 
on  Open 
Boundary7 


Prepare 
Node- Element 
Connectivity 


Prepare  Data 
to  Start  PT 
Computation 


Compute  PT  Using  RK  in 
the  Active  Element 


EBE-based  PT  w 


Update  Active  Element 
Information  tor 
Successive  PT 
Compulation 


Yes 


Particle  Exits: 
PT  Ends 


PT  Ends 


Figure  1.  Element-by-Element  (EBE)  particle  tracking  diagram. 


As  shown  in  Figure  1,  PT123  reads  domain  geometry,  velocity,  and  neces¬ 
sary  information  for  particle  tracking.  It  uses  the  information  of  domain 
geometry  to  prepare  node-element  connectivity,  where  the  elements  con¬ 
necting  at  each  global  node  are  identified  and  stored.  To  track  a  particle, 
PT123  first  locates  the  element  where  the  particle  has  entered.  This 
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element  is  called  the  active  element.  Then  it  conducts  PT  computation 
within  this  active  element  using  the  designated  RK  scheme  (if  the  particle 
is  on  an  interface  between  elements,  all  elements  owning  this  particle  are 
potential  active  elements  and  will  be  examined  one  by  one  until  a  suc¬ 
cessful  PT  computation  is  performed). 

PT123  uses  the  user-specified  initial  time  step  size  for  the  first  PT  com¬ 
putation  within  the  active  element.  It  can  also  compute  a  Courant  number- 
based  time  step  size  for  the  first  PT  computation  in  the  active  element  by 
using  the  CR  parameter  value  specified  in  the  PT  Specific  File  (Appen¬ 
dix  B).  This  Courant  number-based  time  step  size  depends  on  both  the  size 
and  velocities  associated  with  the  active  element.  Suppose  L  represents  the 
characteristic  length  of  the  active  element  which  has  N  element  nodes; 
Vi(fi)  and  Vi(f2)  represent  the  velocity  associated  with  times  ti  and  f2, 
respectively  at  the  z-th  element  node,  where  U  and  f2  are  two  consecutive 
time  steps;  and  the  L  is  the  length  of  a  l-D  element,  the  square  root  of  the 
area  of  a  2-D  element,  and  the  cube  root  of  the  volume  of  a  3-D  element. 
Then  the  Courant  number-based  time  step  size  is  computed  using 
Equation  11. 


At 


Courant 


CRx  — 
V 

avg 


(li) 


where: 

CR  =  name  of  the  Courant  number  parameter 
V  =  average  element  velocity,  which  can  be  computed  using 


N  N 

EIV.M+EIV.M 


v 


2=1 


avg 


2 -N 


(12) 


PT123  prevents  A tCourant  from  becoming  too  large  by  restricting  Vavg  to  a 
minimum  value  of  to-10. 

Reduction  of  time  step  size  at  the  element  boundary 

If  the  time  step  size  used  is  too  large  such  that  the  particle  will  go  outside 
the  active  element,  the  time  step  size  will  be  reduced  so  that  the  particle 
would  reach  the  boundary  of  the  active  element.  This  reduction  of  time 
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step  size  is  enforced  in  the  EBE  tracking  even  if  a  non-adaptive  RK  scheme 
is  employed.  Figure  2  demonstrates  how  this  reduction  of  tracking  time 
step  size  is  achieved  in  PT123. 


Figure  2.  Plot  to  demonstrate  how  PT  time  step 
size  is  reduced  when  the  end  location 
is  outside  of  the  active  element. 


In  Figure  2,  points  S  and  E  represent  the  start  and  the  end  locations  of  a 
PT  computation  associated  with  the  active  element  M,  and  point  I  is  the 
intercept  of  segment  SE  and  the  element  boundary.  If  the  PT  time  step  size 
from  point  S  to  point  E  is  A t0id,  then  the  new  PT  time  step  size  to  prevent 
the  particle  from  going  outside  the  element  is  estimated  using  the 
following  equation: 


A  t 


new 


At  i 

LXLold  1 

\SE 


(13) 


where: 

At mw  =  new  PT  time  step 

lsl  =  distance  between  points  S  and  I 
ZSE  =  distance  between  points  S  and  E. 

It  is  noted  that  this  linear  reduction  of  PT  time  step  size  (i.e.,  Equation  13) 
will  get  the  particle  to  point  E  if  RKi  is  used.  When  higher-order  RK 
schemes  are  used,  it  would  require  an  iteration  process  to  stop  the  particle 
on  the  element  boundary.  To  help  stop  the  particle  on  the  boundary  of  the 
active  element  without  spending  too  computation  in  the  iteration  process, 
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PT123  employs  a  narrow  buffer  zone  surrounding  the  element  boundary. 
When  the  end  location  of  a  PT  computation  is  within  the  buffer  zone,  it  is 
considered  on  the  element  boundary  and  this  end  location  is  adjusted  to 
assure  that  the  particle  is  right  on  the  element  boundary  for  the  subse¬ 
quent  PT.  This  narrow  buffer  zone  is  small  when  compared  to  the  size  of 
the  active  element.  This  buffer  zone  is  defined  using  the  parameter 
DN_SAFE  specified  in  the  PT123  Super  File  (Appendix  B)  and  the  absolute 
error  tolerance,  i.e.,  ATOL,  mentioned  in  Equation  7.  For  example,  if  a 
tracked  particle  is  to  exit  the  active  element  M  via  the  boundary  side  2-3 
as  shown  in  Figure  2,  the  interpolation  factor  associated  with  element 
node  1,  say  DNi,  at  the  exit  location  will  be  zero.  During  PT  computation, 
when  DNi  at  the  end  location  is  computed  between  (o  -  DN_SAFEi )  and 
(o  +  DN_SAFEi),  PT123  considers  this  particle  reaching  the  element 
boundary.  And  then  PT123  sets  DNi  to  zero  and  adjusts  the  interpolation 
factors  associated  with  nodes  2  and  3,  i.e.,  DN2  and  DN3,  accordingly  to 
move  the  particle  slightly  onto  the  element  boundary.  The  buffer  zone 
parameter  DN_SAFEi  is  defined  as  follows. 

1 0 x  ATOT 

DN  _  SAFE1  =  DN  _  SAFE  +  ( 14) 

L 

where  L  is  the  characteristic  length  of  the  active  element  M,  as  mentioned 
before  in  Equation  11. 

Both  the  end  location  and  the  available  tracking  time  are  examined  after 
each  successful  PT  segment  is  computed.  If  the  end  location  is  still  within 
the  active  element  and  the  available  tracking  time  is  not  zero,  successive 
PT  computations  are  conducted  until  either  the  tracking  time  is  completely 
consumed  or  the  particle  reaches  the  boundary  of  the  active  element. 

The  cumulative  tracking  of  a  particle  is  considered  complete  when  either 
the  available  tracking  time  becomes  zero  or  the  particle  exits  from  an  open 
boundary  before  tracking  time  is  consumed  completely.  An  open  boundary 
is  a  boundary  through  which  particles  are  permitted  to  enter  or  leave  the 
domain  of  interest.  When  the  particle  reaches  the  boundary  of  the  active 
element  that  is  not  an  open  boundary  and  the  available  tracking  time  is 
not  zero,  tracking  will  continue  on.  In  this  case,  all  active  element  candi¬ 
dates  are  tested  one  by  one  until  a  successful  tracking  is  conducted  as 
mentioned  before.  This  thus  forms  the  EBE-based  tracking  as  highlighted 
with  the  yellow  shade  in  Figure  1. 
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2.5  Non-Element-by-Element  (NEBE)  tracking 

PT123  can  also  conduct  PT  on  a  NEBE  basis,  where  each  tracking  segment 
computed  may  be  across  elements.  The  NEBE-based  tracking  in  PT123 
uses  the  designated  RK  scheme  to  compute  the  estimated  end  location, 
i.e.,  xn+i,  and  a  ray  search  algorithm  to  locate  the  element  that  includes 
xn+i.  This  process  continues  until  PT  calculations  are  completed  over  the 
entire  computational  domain.  The  PT  time  step  size  will  be  reduced  when¬ 
ever  the  particle  encounters  the  domain  boundary.  Figure  3  depicts  the 
computational  flow  chart  using  the  NEBE-based  PT. 


Figure  3.  Non-Element-by-Element  (NEBE)  particle  tracking  diagram. 


As  shown  on  the  upper  left  of  the  green-shade  area  in  Figure  3,  the  PT 
time  step  size  is  to  be  reduced  if  the  particle  would  go  outside  of  the 
domain  of  interest,  where  no  velocity  field  is  available.  The  computation  of 
the  reduced  time  step  is  given  as  follows  when  various  RK  schemes  are 
used. 
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2.5.1  When  using  RK1  for  NEBE-based  PT: 

The  estimated  end  location  can  be  computed  using 

x„+i=x„+AVv(x„>0  +  0(Afo2)  (!5) 

If  the  estimated  end  location,  i.e.,  xn+i,  is  outside  of  the  domain  of  interest, 
and  the  line  segment  connecting  xn  and  xn+i  intersects  with  the  domain 
boundary  at  xb,  then  the  reduced  PT  time  step  size  is  computed  as 

At1  =  Af0  •  lXfi~X"l  (16) 

x„^-x„ 


where: 

A t±  =  new  PT  time  step  after  reduction. 

2.5.2  When  using  RK2  for  NEBE-based  PT: 

The  estimated  end  location  can  be  computed  using 

k1=Af-V(xn,tn) 

k2  =  AfV(xn+|k1,t„+|At)  (17) 

X„+l=Xn+k2+°(Af3) 

There  are  two  situations  where  the  time  step  size  reduction  may  be  needed 
during  the  RK2  computation.  The  first  instance  occurs  when  the  computed 
xn  +  (ki/2)  is  outside  of  the  domain  of  interest,  and  the  line  segment  con¬ 
necting  x„  and  xn  +  (ki/2)  intersects  with  the  domain  boundary  at  xbi.  In 
this  case,  the  PT  time  step  size  is  reduced  using 

A7,=At0lX”~^”l  (18) 

When  xn  +  (ki/2)  is  within  the  domain,  k2  can  be  estimated.  In  the  second 
case,  if  the  estimated  end  location,  i.e.,  xn+i,  is  outside  of  the  domain  of 
interest,  and  the  line  segment  connecting  xn  and  xn+i  intersects  with  the 
domain  boundary  at  xs2,  then  the  PT  time  step  size  is  reduced  using 
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At±  =  At0  -  ,  ^  X"  (19) 

|Xn+l-X„| 

2.5.3  When  using  RK4  for  NEBE-based  PT: 

The  estimated  end  location  can  be  computed  using 
k1=A  f-V(xn,fn) 

k2  =  At- V(xn  +|k15tn  +|a  t) 

k3=At-V(xn+|k2,tn+|At)  (20) 

k4  =  At- V(xn  +k3,tn  +  At) 

Xn+1  =  Xn  +  +^k4  +0(At5) 

6  3  3  6 

There  are  four  possible  situations  for  the  reduction  of  PT  time  step  size 
during  the  RK4  computation.  These  four  conditions  are  described  next. 

Situations  1-2:  the  computed  xn  +  (k,/2)  (z  =  1,  2)  is  outside  of  the 
domain  of  interest,  and  the  line  segment  connecting  xn  and  xn  +  (k,/2) 
intersects  with  the  domain  boundary  at  xbu  In  this  case,  the  PT  time  step 
size  is  reduced  using 


A  K 


A  tn 


|Xj?z  Xr 

Ik../2| 


(21) 


Situation  3:  the  computed  xn  +  k3  is  outside  of  the  domain  of  interest, 
and  the  line  segment  connecting  xn  and  xn  +  k3  intersects  with  the  domain 
boundary  at  xs3.  In  this  case,  the  PT  time  step  size  is  reduced  using 

At.  =Af0-t*AlZXJ  (22) 

k. 


Situation  4:  the  estimated  end  location,  i.e.,  xn+i,  is  outside  of  the 
domain  of  interest,  and  the  line  segment  connecting  xn  and  xn+i  intersects 
with  the  domain  boundary  at  xs4.  Here  the  PT  time  step  size  is  reduced 
using 
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At±  =  At0  •  }Xb4  X"  (23) 

Ixn+i-x„| 

2.5.4  When  using  RK45  for  NEBE-based  PT: 

The  estimated  end  location  can  be  computed  using 

k,  =  At- V(xn,tn) 

k2  =  At  ■  V(xn  +  621k1 ,  tn  +  a2At) 

k3  =  At  •  V(xn  +  b31  k,  +  h32k2 ,  tn  +  a3At) 

k4  =  At  ■  V (x  +  5  k4  +  b42k2  +  543k3  ,t  +  a4Af) 

(24) 

k5  =  At-  V(xn  +651k4  +h52k2  +b53k3  +h54k4,fn  +a5Af) 
k6  =  At-  V(xn  +b61ls.1  +  b62k2+b63k3+b64k4  +b65k5,tn+a6At) 
xni  =  x„  +  c1k1  +  c2k2  +  c3k3  +  c4k4  +  c5k5  +  c6k6  +  0(Af6) 
xLi  =  xn  +  ciki  +  c2k2  +  c3k3  +  c4k4  +  c*k5  +  c*k6  +  0(A  f 5 ) 

There  are  seven  possible  situations  for  the  reduction  of  PT  time  step  size 
during  the  RK45  computation.  These  seven  situations  are  described  as 
follows. 


Situations  1-5:  the  computed  xn  +  y~]bfa-k,-  ( k  =j+ 1,  i  =  1,  2,  3,  4,  5)  is 

j= 1 

outside  of  the  domain  of  interest,  and  the  line  segment  connecting  xn  and 

i 

xn  +  y^hb-k.  intersects  with  the  domain  boundary  at  xbu  In  this  case,  the 

j= 1 

PT  time  step  size  is  reduced  using 


At1=At0- 


j= 1 


(25) 


Situation  6:  the  computed  xn+i  is  outside  of  the  domain  of  interest,  and 
the  line  segment  connecting  xn  and  xn+i  intersects  with  the  domain  boun¬ 
dary  at  xb6-  Here  the  PT  time  step  size  is  reduced  using 


At±  —  At0  ■ 


|X,B6  Xn 

|Xn+l  Xr 


(26) 
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Situation  7:  the  computed  x^+1  is  outside  of  the  domain  of  interest,  and 
the  line  segment  connecting  xn  and  xjn+1  intersects  with  the  domain  boun¬ 
dary  at  xb7.  In  this  case,  the  PT  time  step  size  is  reduced  using 


At,  =  At0 


XB7~Xr. 


X„+1-Xr 


(27) 


When  the  reduction  of  PT  time  step  size  happens,  A t0  is  updated  with  A t± 

to  compute  xn+i  using  the  specified  RK  scheme.  Also,  the  same  examina¬ 
tion  of  whether  the  PT  time  step  size  needs  to  be  further  reduced  is  con¬ 
ducted.  This  process  continues  until  the  computed  xn+i  is  within  the 
domain. 

2.6  Tracking  along  a  closed  boundary 

An  open  boundary  is  a  boundary  through  which  a  particle  can  enter  or  exit 
the  domain  of  interest.  A  boundary  is  a  closed  boundary  if  it  is  not  an  open 
boundary.  Conceptually,  the  flow  velocity  associated  with  a  closed  boun¬ 
dary  is  parallel  or  tangential  to  the  boundary,  i.e.,  zero  normal  velocity  at 
the  closed  boundary.  However,  both  mesh  resolution  and  numerical  error 
can  contribute  to  non-tangent  flow  velocity  at  the  closed  boundary.  This  is 
common  in  the  simulations  of  real-world  problems.  As  a  result,  the  com¬ 
puted  PT  results  can  become  misleading  if  the  PT  computation  does  not 
proceed  when  the  tracked  particle  reaches  a  closed  boundary. 

PT123  uses  the  projected  velocity  on  the  closed  boundary  to  continue  PT 
computation  until  the  tracked  particle  reaches  an  open  boundary  or  when 
the  tracking  time  is  completely  consumed.  The  computation  of  projected 
velocity  on  the  closed  boundary  is  described  next. 

2.6.1  Velocity  projection  onto  a  2-D  boundary  edge 

As  shown  in  Figure  4,  the  computed  velocity  at  node  1  is  Vi,  which  is  non¬ 
parallel  to  the  closed  boundary  edge  between  nodes  1  and  2.  The  geometric 
quantities  associated  with  edge  1-2  for  PT123  computations  include: 
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VI  (Vxl,  Vyl) 


2  (x2,  y2) 


1  (xl.yl) 


Figure  4.  Projection  of  velocity  onto  a  2-D  boundary  segment. 


l.  The  length  of  edge  1-2  (l12): 


(28) 


where: 

Ax  =  X2  -  xi 
Ay  =  yi  -  yi 

(xi,  yi)  =  coordinates  of  node  1 
(x2,  y  2)  =  coordinates  of  node  2. 

2.  The  unit  vector  parallel  to  edge  1-2,  u,  is 


(29) 


3.  The  projected  magnitude  of  Vi  onto  edge  1-2  can  be  computed  as 


(30) 


Vip  =  Vi  u 


where: 

Vip  =  projected  velocity  of  Vi  onto  edge  1-2. 

Then  the  projected  velocity  of  Vi  onto  edge  1-2,  i.e.,  Vip,  is  computed 


as 
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Vip  =  |Vip|u  =  (Vi'u)u 


Tr  .  Ax  Tr  .  Ay 
Vxl - hVyl-  — 

l12  l. 


12 


Ax  Ay 

~’T~ 
l12  l12 


(31) 


where: 

Vxl,  Vy  1  =  x-  and  y-components  of  Vi. 

2.6.2  Velocity  projection  onto  a  3-D  boundary  face 

Figure  5  shows  the  geometric  relationship  of  the  velocity  at  node  1,  i.e.,  Vi, 
and  a  3-D  triangular  boundary  face  1-2-3.  The  equation  describing  the 
plan  containing  face  1-2-3  can  be  represented  by  ax  +  by  +  cz  +  d  =  o, 
where  the  normal  vector  of  the  plane  is  (a,  b,  c). 


The  normal  vector  of  the  plane,  n,  can  be  computed  using  Equation  32  as 
n  =  (a,A,c) ^AL12,AL13  =  (Ax12,Ay12,Az12)x(  x13  y13  z13)  (32) 


CM 

t-H 

<1 

AZi2 

AFl3 

^13 

AZi2 

AXi2 

^13 

^13 

where: 
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C 


AX^ 

ky±2 


A z 
Ax 


12 

13 


Ayi3 


(xl,  yl,  zl) 
(x2,  y2,  z2) 
(x3,  y3,  z3) 


Ax12  Ay. 

Ax13  Ay, 

x2  — 

xl 

y2- 

yi 

z2  - 

zl 

x3  - 

xl 

y3- 

yi 

z3  - 

zl 

coordinates  of  node  l 
coordinates  of  node  2 
coordinates  of  node  3. 


The  unit  normal  vector  is  calculated  as 

a  b  c 

u=  i — r»i — r»i — i 

n  n  n 


(33) 


where: 


u  =  the  unit  normal  velocity  of  the  plan  containing  face  1-2-3 

|n|  =  yja2  +  b2  +c2  . 

The  projected  velocity  of  Vi  parallel  to  the  unit  normal  velocity  u  is 


|Vin  u  =  (Vi-u)u  = 

\ 

Vxl-  a  +  Vyl  Q  +  Vzl-  Q 

\ 

a  a  a 

n  n  n 

/ 

n  n  n 

/ 

(34) 


where: 

Vxl,  Vyl,  Vzl  =  x- ,  y-,  and  z-components  of  Vi 

The  projected  velocity  of  Vi  onto  Face  1-2-3  is 
Vip  =  Vl-Vin  =  Vi-|Vin|u 

=  (Vxl,  Vyl,  Vzl)  - 


\ 

ci  b  c 

\ 

a  b  c 

+  Vyl-  +  Vyl- 

n  n  n 

/ 

n  n  n 

/ 

(35) 
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It  is  important  to  note  that  each  boundary  element  should  have  only  one 
closed  boundary  edge/face  when  the  unstructured  mesh  is  constructed. 
Otherwise,  the  velocity  at  nodes  associated  with  two  closed  boundary 
edges/faces  will  not  be  projected  properly. 
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3  Test  Examples 

This  chapter  presents  six  PT  examples  using  PT123.  The  first  example  was 
designed  to  compare  the  use  of  node-based  and  element-based  velocities 
for  a  l-D  heterogeneous  problem.  Examples  2-4  were  employed  for  veri¬ 
fication  and  comparison  among  different  tracking  schemes.  The  last  two 
examples  demonstrate  the  application  of  PT123  in  real-world  problems. 
Because  the  first  four  examples  were  designed  mainly  to  examine  PTi23’s 
numerical  techniques,  any  time  and  length  units  can  be  employed.  For 
convenience  in  discussion,  we  simply  used  meter  (m)  and  second  (sec)  as 
the  length  and  time  units,  respectively,  for  these  four  examples. 

3.1  Example  1:  l-D  steady  non-uniform  velocity  field 

In  this  example,  the  l-D  domain  was  composed  of  41  nodes  and  40  linear 
elements.  Every  element  had  the  same  length  of  10  m.  The  domain 
included  four  types  of  material:  material  type  1  for  elements  1-5,  8-11, 
15-17,  25-34;  material  type  2  for  elements  6-7  and  35-40;  material  type 
3  for  elements  12-14  and  22-24;  and  material  type  4  for  elements  18-21. 
The  materials  have  distinct  levels  of  permeability  for  water  flow.  By 
employing  a  steady  flow  throughout  the  entire  domain  (from  left  to  right), 
different  velocities  appeared  in  distinguishable  parts  of  the  domain 
(Figure  6)  due  to  heterogeneity.  To  highlight  the  heterogeneity,  each  linear 
element  in  Figure  6  was  depicted  with  a  rectangle  filled  with  a  color  repre¬ 
senting  the  material  type.  In  Figure  6,  both  element-based  and  node-based 
velocities  are  provided.  The  element-based  velocity  was  computed  based 
on  the  material  type  associated  with  each  element,  where  higher  velocities 
were  in  the  elements  associated  with  material  types  that  are  more 
permeable. 


Q, 


Element-based  Velocity: 

2  0.2  2  0.02  2  20  0.02 
h - - >l<  >K 


0.2 


Q 


2  1 0.2  \  2  I  0.02  \  2  \  20  j  0.02  \  2  j  0.2 

1.1  1.1  1.01  1.01  11  11.01  1.01  1.1 
Node-based  Velocity: 

□  Material  Type  1  P  Material  Type  2 
P  Material  Type  3  P  Material  Type  4 


Figure  6.  Element-  and  node-based  velocity  variation  for  Example  1. 
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On  the  other  hand,  the  node-based  velocity  was  determined  from  element 
connectivity  with  proportional  contribution.  As  a  result,  a  node  at  the 
interface  of  two  material  types  has  an  average  velocity. 

To  test  PT123, 10  particles  were  populated  at  the  center  of  elements  1,  3,  5, 
7,  9, 11, 13, 15, 17,  and  19,  respectively.  A  PT  computation  over  a  time 
period  of  too  sec  was  conducted  using  both  element-based  and  node- 
based  velocities,  and  both  element-by-element  (EBE)  and  non-element- 
by-element  (NEBE)  tracking  methods.  Although  most  existing  flow  models 
compute  node-based  velocity,  the  velocity  field  specification  type  does  not 
represent  the  velocity  change  on  the  interface  of  two  material  types  accu¬ 
rately.  Instead,  element-based  velocity  can  represent  the  true  velocity  field 
correctly.  For  this  simple  l-D  example,  the  analytical  solution  of  particle 
tracking  can  be  obtained  using  the  element-based  velocity  according  to 
Darcy’s  law  (Jury  et  al.  1991),  where  the  time  duration  expended  by  a 
particle  to  pass  through  an  element,  A tM ,  is  equal  to 

(36) 


where: 

Lm  =  the  length  of  line  element  M 

VM  =  the  element-based  velocity  of  element  M. 


Table  2  lists  the  mean  absolute  error  (MAE)  associated  with  each  particle 
when  the  PT123  results  are  compared  with  the  analytical  solution  every 
second  during  the  computation.  The  MAE  for  a  particle  is  defined  as 


MAE 


computed 


^  j  analytical 


N 


(37) 


where: 


N  =  number  of  comparisons  (it  is  too  here  for  the  comparison  is 
made  every  second) 

xcomputed  (£_)  _  compUteci  particle  location  at  the  time  associated  with  the 
i-th  comparison  (it  is  at  time  =  i  sec) 
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x'maiyticai  )  _  anaiytical  particle  location  at  the  time  associated  with  the 
i-th  comparison  (it  is  at  time  =  i  sec). 


Table  2.  Mean  absolute  errors  for  Example  1. 


MAE 

Particle  ID 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

EBE_l_e_0.1 

<  IE-6 

EBE_4_e_0.1 

<  IE-6 

EBE_l_e_l 

<  IE-6 

EBE_4_e_l 

<  IE-6 

NEBE_l_e_0.1 

0.139 

0.156 

0.174 

0.045 

0.172 

0.192 

0 

0.204 

0.223 

0.979 

NEBE_4_e_0.1 

0.023 

0.026 

0.029 

0.043 

0.028 

0.032 

0 

0.601 

0.670 

0.652 

NEBE_l_e_l 

0.695 

0.784 

0.891 

0.454 

0.862 

0.960 

~0 

1.021 

1.119 

14.68 

NEBE_4_e_l 

0.463 

0.522 

0.585 

0.435 

0.575 

0.640 

~0 

0.680 

0.746 

4.896 

EBE_l_n_0.1 

16.75 

21.91 

27.15 

13.19 

6.956 

7.818 

0.5 

9.707 

10.19 

9.609 

EBE_4_n_0.1 

16.73 

21.89 

27.14 

13.21 

6.948 

7.811 

0.5 

9.756 

10.20 

9.597 

EBE_l_n_l 

16.93 

22.10 

27.29 

13.04 

7.024 

7.886 

0.5 

9.390 

10.18 

9.757 

EBE_2_n_l 

16.69 

21.85 

27.10 

13.19 

6.946 

7.808 

0.5 

9.576 

10.11 

9.510 

EBE_4_n_l 

16.73 

21.89 

27.14 

13.21 

6.948 

7.810 

0.5 

9.752 

10.20 

9.587 

NEBE_l_n_0.1 

16.57 

21.94 

27.15 

13.19 

6.956 

7.818 

0.5 

9.706 

10.19 

9.609 

NEBE_4_n_0.1 

16.73 

21.89 

27.14 

13.21 

6.948 

7.811 

0.5 

9.756 

10.20 

9.597 

NEBE_l_n_l 

16.93 

22.10 

27.30 

13.04 

7.025 

7.888 

0.5 

9.314 

10.14 

9.792 

NEBE_2_n_l 

16.70 

21.85 

27.10 

13.20 

6.946 

7.808 

0.5 

9.596 

10.09 

9.482 

NEBE_4_n_l 

16.74 

21.89 

27.14 

13.21 

6.948 

7.811 

0.5 

9.764 

10.19 

9.593 

In  Table  2,  the  PT123  simulations  are  listed  with  names  indicating  the  four 
test  variables  specified:  “EBE”  and  “NEBE”  identify  the  tracking  method; 
the  first  number  indicates  the  RK  scheme;  “e”  and  “n”  denote  the  use  of 
element-based  and  node-based  velocity,  respectively;  and  the  second 
number  is  the  time  step  size  for  integration.  For  example,  simulation 
“EBE_i_e_o.i”  used  the  element-by-element  tracking  method,  the 
ist-order  RK  scheme,  an  element-based  velocity  field,  and  a  time  step  size 
of  0.1. 
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Figures  7-10  compare  the  PT123  results  using  various  RK  schemes  with 
the  analytical  solution  associated  with  particles  5  and  10,  where  the 
element-based  velocity  was  implemented  in  Figures  7  and  8,  and  the  node- 
based  velocity  was  applied  in  Figures  9  and  10.  The  legends  of  Figures  7- 
10  refer  to  the  simulations  using  the  same  naming  system  explained  above. 


Particle  5  Trajectory  (Using  Element-Based  Velocity) 


Particle  10  Trajectory  (Using  Element-Based  Velocity) 
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Figure  7.  Comparison  of  the  tracking  paths  of  particles  5  and  10  using  element-based  velocity 

and  various  EBE  tracking  schemes. 


Particle  5  Trajectory  (Using  Element-Based  Velocity) 


Particle  10  Trajectory  (Using  Element-Based  Velocity) 


Figure  8.  Comparison  of  the  tracking  paths  of  particles  5  and  10  using  element-based  velocity 

and  various  NEBE  tracking  schemes. 


Particle  5  Trajectory  (Using  Node-Based  Velocity) 


Particle  10  Trajectory  (Using  Node-Based  Velocity) 
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Figure  9.  Comparison  of  the  tracking  paths  of  particles  5  and  10  using  node-based  velocity 

and  various  EBE  tracking  schemes. 
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Particle  5  Trajectory  (Using  Node-Based  Velocity) 


Particle  10  Trajectory  (Using  Node-Based  Velocity) 
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Figure  10.  Comparison  of  the  tracking  paths  of  particles  5  and  10  using  node-based  velocity 

and  various  NEBE  tracking  schemes. 


Figure  7  shows  that  all  four  computational  results  from  EBE-based  track¬ 
ing  can  be  considered  to  match  the  analytical  solutions  of  particles  5  and 
10  perfectly.  This  illustrates  the  advantage  of  using  EBE-based  tracking 
when  the  element-based  velocity  is  employed.  Figure  8  and  Table  2,  on  the 
other  hand,  show  that  the  larger  the  PT  time  step  size  is,  the  greater  error 
NEBE-based  tracking  produces.  The  error  associated  with  NEBE-based 
tracking  grows  with  the  discontinuity  of  velocity  from  one  material  type  to 
another  (e.g.,  from  material  type  4  to  material  type  3).  It  thus  suggests  that 
NEBE  cannot  provide  accurate  tracking  results  when  element-based 
velocity  field  is  used  for  domains  with  heterogeneity. 

Figures  9  and  10  show  that  both  EBE-  and  NEBE-based  tracking  methods 
produce  similar  results  for  particles  5  and  10.  They  also  reveal  the  tracking 
error  resulting  from  implementing  node-based  velocity  when  compared 
with  the  analytical  solution.  As  shown  in  Table  2,  when  the  node-based 
velocity  was  used,  i.e.,  the  bottom  10  test  cases,  the  tracking  error  remains 
significant  no  matter  which  RK  scheme  or  time  step  size  were  specified.  It 
is  noted  that  the  error  associated  node-based  velocity  can  be  reduced  by 
including  small  elements  around  the  material  interface,  but  cannot  be 
completely  removed. 

3.2  Example  2:  2-D  steady  rotational  velocity 

In  this  example,  PT  was  computed  in  a  2-D  square  domain,  ranging  from 
-2,000  m  to  2,000  m  in  both  the  x-  and  the  y-directions.  The  domain  was 
discretized  using  both  quadrilateral  and  triangular  elements  (Figure  11). 
The  mesh  was  composed  of  81  nodes  and  96  elements.  A  steady  velocity 
vector  was  given  at  each  node  using  Equations  38  and  39  (Figure  11). 
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Vx(xyy)j=  0.0Q-2 
Vy(xjy)x=  -O.J0O2 


G  [-2000,2000],  G  [—2000,2000]  (38) 

G  [-£000,2000],  G  [-2000,2000]  (39) 


Figure  11.  The  mesh  and  steady  rotational 
velocity  field  of  Example  2. 


Four  hundred  particles  (point  ID’s  o  -  399)  were  placed  uniformly  on  a 
circle  centered  at  (0,1000)  with  a  radius  of  700  m  (Figure  12).  In  Fig¬ 
ure  12,  each  particle  is  connected  to  its  adjacent  two  particles  with  line 
segments  to  form  a  400-edge  polygon.  With  the  given  steady  velocity  field, 
each  particle  should  move  clockwise  around  the  center,  i.e.,  (0,0),  and 
return  to  its  initial  position  at  times  that  are  multiples  of  1,000  sec  if 
tracking  is  accurate  (Pokrajac  and  Lazic  2002;  Cheng  et  al.  1996). 


Figure  12.  Four  hundred  particles  at  time  =  0  sec 
in  Example  2. 
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The  results  of  PT  simulations  using  different  RK  schemes  and  time  step 
sizes  were  compared  as  shown  in  Table  3  and  Figure  13. 


Table  3.  Example  2  efficiency  comparison. 


Method 

TSS1 

Ht, a 

nt,t 

ei  (m) 

eoo  (m) 

EBE_RK12 

0.002 

200,005,167 

200,020,212 

0.4045E-01 

0.6715E-01 

EBE_RK12 

0.01 

40,005,109 

40,020,286 

0.2022E+00 

0.3355E+00 

EBE_RK12 

0.1 

4,005,182 

4,020,295 

0.2023E+01 

0.3356E+01 

EBE_RK12 

0.5 

805,422 

820,929 

0.1013E+02 

0.1681E+02 

NEBE_RK12 

0.002 

200,000,00 

200,000,000 

0.4044E-01 

0.6711E-01 

NEBE_RK12 

0.01 

40,000,000 

40,000,000 

0.2022E+00 

0.3356E+00 

NEBE_RK12 

0.1 

4,000,000 

4,000,000 

0.2024E+01 

0.3358E+01 

NEBE_RK12 

0.5 

800,000 

800,000 

0.1016E+02 

0.1686E+02 

EBE_RK22 

0.1 

4,006,687 

4,025,247 

0.9681E-03 

0.1934E-02 

EBE_RK22 

1 

409,283 

432,963 

0.4250E-01 

0.7050E-01 

EBE_RK22 

10 

54,861 

85,843 

0.3906E+01 

0.6478E+01 

EBE_RK22 

100 

26,604 

66,423 

0.1026E+03 

0.1958E+03 

NEBE_RK22 

0.1 

4,000,000 

4,000,000 

0.9284E-03 

0.1540E-02 

NEBE_RK22 

1 

400,000 

400,000 

0.4285E-01 

0.7111E-01 

NEBE_RK22 

10 

40,000 

40,000 

0.4235E+01 

0.7029E+01 

NEBE_RK22 

100 

4,000 

4,000 

0.4218E+03 

0.7317E+03 

EBE_RK42 

10 

54,320 

84,581 

0.8066E-04 

0.4005  E-03 

EBE_RK42 

25 

33,597 

67,792 

0.2345E-01 

0.4160E-01 

EBE_RK42 

100 

26,672 

65,833 

0.1101E+01 

0.3948E+01 

EBE_RK42 

500 

26,288 

66,052 

0.4904E+01 

0.4815E+02 

NEBE_RK42 

10 

40,000 

40,000 

0.3328E-03 

0.5523E-03 

NEBE_RK42 

25 

16,000 

16,000 

0.3214E-01 

0.5334E-01 

NEBE_RK42 

100 

4,000 

4,000 

0.8311E+01 

0.1379E+02 

NEBE_RK42 

500 

1,600 

1,600 

0.2845E+03 

0.4721E+03 

EBE_RK452 

100 

1,347,292 

1,456,129 

0.4188E-03 

0.1062E-02 

EBE_RK452 

500 

1,347,351 

1,456,707 

0.4191E-03 

0.1062E-02 

EBE_RK452 

Cr3  =  1 

1,347,359 

1,456,676 

0.4190E-03 

0.1062E-02 

NEBE_RK452 

100 

1,366,565 

1,369,315 

0.4088E-03 

0.6786E-03 

NEBE_RK452 

500 

1,366,568 

1,369,354 

0.4089  E-03 

0.6786E-03 

iTSS  =  time  step  size. 

2  ATOL  =  10-7,  DN_SAFE  =  10  7. 

3  Cr  =  Courant  number. 
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Figure  13.  Mean  absolute  error  versus  computational  WORK  (top) 
and  time  step  size  (bottom)  for  various  RK  schemes  using  EBE 
and  NEBE  tracking  for  Example  2. 


Table  3  provides  a  comparison  of  the  efficiency  for  various  RK  schemes 
employed  in  both  the  EBE-  and  the  NEBE-based  PT.  Here,  nt,t  is  the  total 
number  of  tracking  time  steps  attempted,  nt,a  is  the  total  number  of  steps 
accepted,  and  “Cr  =  CR”  refers  to  choosing  the  initial  time  step  for  tracking 
over  an  element  using  a  local  target  Courant  number  of  CR,  as  defined  in 
Equation  11.  When  an  adaptive  RK  scheme  is  used  and  the  estimated  error 
is  greater  than  the  prescribed  error  tolerance,  the  attempted  PT  is  not  a 
successful  segment.  The  attempted  PT  becomes  a  successful  segment  when 
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time  step  size  is  reduced  to  a  degree  such  that  the  prescribed  accuracy  is 
satisfied.  Another  example  of  an  unsuccessful  tracking  segment  is  when 
the  EBE-based  PT  is  employed  and  a  particle  would  go  outside  the  active 
element,  time  step  size  needs  to  be  reduced  so  that  the  particle  would 
reach  the  boundary  of  the  active  element.  This  is  why  when  EBE-based  PT 
or  RK45  is  used,  nt,t  can  be  greater  than  nt,a.  Two  types  of  error  indica¬ 
tion  are  also  included  in  Table  3  for  accuracy  comparison,  where  errors 
were  measured  using  the  analytical  solution  (Cheng  et  al.  1996)  at  times 
tn  =  i,ooon  sec,  n  =  1, ...  10.  These  two  error  indicators  are  defined  as 
follows. 


ei 


1  1 
~10X  NPTX 


E 


i=l,NPT 
n= 1,10 


(40) 


e  =  max 

00  i=l,NPT 


xr(fn)-xr*“(fn) 


(41) 


where: 

NPT  =  number  of  particles. 

To  compare  the  computational  effort  for  each  test  case  more  closely,  we 
define  a  number  WORK  as  follows. 

WORK  =  nt  t  x  Nstage  (42) 

where  Nstageis  equal  to  1  for  RKi,  2  for  RK2,  4  for  RK4,  and  6  for  RK45. 

Figure  13  plots  mean  absolute  error,  i.e.,  s1}  versus  WORK  and  time  step 
size  for  various  RK  schemes  using  EBE  and  NEBE. 

It  must  be  noted  that  the  WORK  associated  with  NEBE-based  tracking,  as 
shown  in  Figure  13,  does  not  include  the  computation  spent  for  ray  trac¬ 
ing.  From  Table  3  and  Figure  13,  the  following  are  observed: 

1.  For  non-adaptive,  NEBE-based  PT,  nt,a  =  nt,t  and  total  tracking  time 
(i.e.,  10,000  sec  in  Example  2)  =  TSS  x  rit,t,  where  TSS  is  time  step  size. 

2.  For  EBE-based  PT,  nt,a  <  nt,t  due  to  the  reduction  of  TSS  when  the  par¬ 
ticle  would  go  outside  of  the  active  element. 
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3.  Given  specified  TSS,  non-adaptive,  EBE-based  PT  yields  smaller  si,  i.e., 
more  accurate,  when  compared  with  the  respective  non-adaptive, 
NEBE-based  PT. 

4.  When  the  higher-order  RK  scheme  is  used,  the  larger  TSS  can  be  used 
to  obtain  accurate  results  when  non-adaptive  PT  is  considered. 

5.  When  the  non-adaptive  tracking  is  considered,  the  EBE-based  PT 
requires  more  RK  computation  (i.e.,  WORK)  than  the  NEBE-based  PT 
due  to  the  reduction  of  TSS  when  the  particle  would  go  outside  of  the 
active  element. 

6.  Given  specified  error  tolerance,  the  tracking  effort  (i.e.,  WORK)  is 
insensitive  to  TSS  when  adaptive  schemes  are  used. 

3.3  Example  3:  2-D  swirl  velocity 

This  example  also  computed  PT  in  a  2-D  square  domain,  though  ranging 
from  o  to  1  m  in  both  the  x-  and  the  y-directions.  The  domain  was  dis¬ 
cretized  using  800  triangular  elements  and  441  nodes.  Transient  nodal 
velocities  were  computed  at  time  =  0.0,  0.5, 1.0, ...,  79.5,  and  80  sec  using 
Equations  43  and  44.  These  equations  describe  a  velocity  field  five  times 
faster  than  the  velocity  field  employed  for  a  linear  advection  problem 
(Problem  E)  in  (Farthing  and  Kees  2009),  where  the  initial  concentration 
disk  underwent  significant  deformation  during  the  transient  simulation. 


Vx(x,gtf)  =  5cm: 


jit 


Vy(x,ffjt)  =  -5mps 


8 
'  jit 


8 


•sin  (2  )-yin2(  ) 

e  [0,1], 

e[o,l] 

(43) 

\ 

esin(2  y)-sin2(  ) 

e  [0,1], 

6  [0,1] 

(44) 

As  shown  in  Figure  14,  the  velocity  vectors  were  counterclockwise  from 
time  =  o  to  4  sec,  clockwise  from  time  =  4  to  12  sec,  and  counterclockwise 
from  time  =  12  to  16  sec  to  complete  a  cycle.  Velocities  were  zero  at 
time  =  4  and  12  sec  when  the  flow  changed  directions. 

For  testing  PT,  400  particles  were  initialized  forming  a  circle,  where  the 
center  was  at  (0.5,0.75)  and  the  radius  was  0.15  m  (Figure  15).  The  hydro- 
dynamic  time  step,  i.e.,  0.5  sec,  was  used  as  the  initial  time  step  size  for 
PT. 
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Figure  14.  Transient  velocity  fields  at  various  times  of  Example  3. 


Figure  15.  Four  huncred  particles  at  time  =  0  sec 
in  Example  3. 


Figures  16  and  17  show  the  PT  results  from  time  =  o  to  8  sec  and  from 
time  =  8  to  16  sec,  respectively.  Using  RK45,  the  initial  circle  re-appeared 
at  time  =  8  and  16  sec  as  accurate  PT  will  yield  even  though  the  swirl-type 
velocity  field  changed  the  relative  locations  of  the  400  particles  drastically 
during  the  tracking  process. 
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Figure  16.  Example  3  particle  tracking  results  at  time  =  0,  2, 4,  6,  and  8  sec  using  RK45. 


Figure  17.  Example  3  particle  tracking  results  at  time  =  8, 10, 12, 14,  and  16  using  RK45. 


As  seen  in  Figure  18,  the  shape  and  location  of  the  initial  circle  was  main¬ 
tained  at  time  =  16,  32, 48,  64,  and  80  sec,  which  correspond  to  the  end  of 
1,  2,  3,  4,  and  5  cycles  of  PT,  respectively,  when  RK45  was  used. 
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Figure  18.  Example  3  particle  tracking  results  at  time  =  0, 16, 32, 48,  64, 

and  80  using  RK45. 


Table  4  compares  the  efficiency  for  various  RK  schemes  employed  in  both 
the  EBE-  and  the  NEBE-based  PT,  where  errors  were  measured  using  the 
analytical  solution  (Farthing  and  Kees  2009)  at  times  tn  =  8n,  n  =  1, ...  10. 
Table  5  compares  the  computation  profile  that  includes  the  four  most 
called  subroutines  associated  with  each  PT  technique  employed  for 
Example  3.  From  Tables  4  and  5,  the  following  observations  are  made: 

1.  For  adaptive  PT,  reducing  ATOL  and/or  DN_SAFE  would  increase 
accuracy. 

2.  For  adaptive  PT,  both  nt,a  and  nt,t  are  sensitive  to  ATOL  when 

TTS  =  0.1,  but  become  less  sensitive  when  Cr  =  0.1  was  employed  to 
determine  TTS. 

3.  For  both  adaptive  and  non-adaptive  PT,  the  computation  involved  in 
ray  tracing  is  significant  (subroutines  EL_INTERSECTi23  and  ES123 
are  called  for  ray  tracing  computation). 

4.  The  ratio  of  ray-tracing  computation  to  RK  computation  increases  with 
TSS  when  non-adaptive  PT  is  considered. 
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Table  4.  Example  3  efficiency  comparison. 


Method 

TSS1 

n  t,a 

n%t 

ei  (m) 

£oo  (m) 

EBE_RK22 

O.OOl 

35,305,792 

43,383,263 

0.7185E-04 

0.4071E-03 

EBE_RK22 

0.002 

21,615,186 

33,870,002 

0.3960E-03 

0.1649E-02 

EBE_RK22 

0.004 

45,019,632 

119,458,841 

0.2776E-02 

0.1451E-01 

NEBE_RK22 

0.001 

32,000,000 

32,000,000 

0.6274E-04 

0.2851E-03 

NEBE_RK22 

0.002 

16,000,000 

16,000,000 

0.5009E-03 

0.2370E-02 

NEBE_RK22 

0.004 

8,000,000 

8,000,000 

0.4037E-02 

0.2102E-01 

EBE_RK42 

0.005 

11,110,354 

20,688,936 

0.2758E-04 

0.2468 E-03 

EBE_RK42 

0.01 

8,720,160 

19,494,241 

0.6537E-04 

0.3515E-02 

EBE_RK42 

0.02 

7,883,456 

19,452,895 

0.2441E-03 

0.7039E-02 

EBE_RK42 

3 

ll 

O 

13,554,443 

22,603,180 

0.2479  E-04 

0.2084E-03 

EBE_RK42 

3 

ll 

o 

cn 

7,771,777 

19,192,392 

0.4727  E-03 

0.4152E-01 

NEBE_RK42 

0.005 

6,400,000 

6,400,000 

0.3895E-04 

0.4226E-03 

NEBE_RK42 

0.01 

3,200,000 

3,200,000 

0.6058E-03 

0.8630E-02 

NEBE_RK42 

0.02 

1,600,000 

1,600,000 

0.9291E-02 

0.6785E-01 

EBE_RK452 

0.1 

8,790,432 

22,187,890 

0.3475  E-04 

0.3120E-03 

EBE_RK454 

0.1 

33,772,766 

61,688,212 

0.2777E-04 

0.2319E-03 

EBE_RK452 

0.5 

8,794,292 

22,201,438 

0.3440E-04 

0.3102E-03 

EBE_RK452 

O 

OJ 

II 

O 

12,193,618 

21,458,297 

0.2590E-04 

0.2221E-03 

EBE_RK455 

O 

01 

II 

O 

8,127,007 

20,353,079 

0.2721E-03 

0.2174E-02 

EBE_RK456 

O 

01 

II 

O 

9,158,432 

22,563,127 

0.9014E-05 

0.7645E-04 

EBE_RK452 

Cr3  =  0.5 

8,780,233 

21,750,165 

0.3465  E-04 

0.2931E-03 

NEBE_RK452 

0.1 

22,199,342 

35,354,014 

0.1364E-04 

0.1189E-03 

NEBE_RK452 

0.5 

22,207,119 

35,441,682 

0.1399E-04 

0.1213E-03 

iTSS  =  time  step  size. 

2  ATOL  =  10-9,  DN_SAFE  =  10  6. 

3  Cr  =  Courant  number. 


4  ATOL  =  10-10,  DN_SAFE  =  106. 

5  ATOL  =  109,  DN_SAFE  =  10  5. 

6  ATOL  =  10-9,  DN_SAFE  =  10  7. 
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Table  5.  Example  3  computation  profile  comparison. 


Method 

TSS1 

Subroutine 

RK124_EBE_PT 

INTRP123 

ELTRAK123 

VEL123 

EBE_RK42 

0.005 

100,552,293 

92,654,981 

68,492,092 

59,174,421 

EBE_RK42 

0.01 

90,221,706 

81,458,098 

58,395,681 

51,233,224 

EBE_RK42 

0.02 

88,027,668 

78,641,648 

55,353,201 

49,121,878 

EBE_RK42 

o 

~cL 

II 

o 

113,517,654 

106,106,754 

80,096,228 

68,311,294 

EBE_RK42 

Cr3  =  0.5 

86,749,218 

77,527,408 

54,676,776 

48,364,434 

Method 

TSS 

Subroutine 

ELJNTERSECT123 

INTRP123 

ES123 

RK4_NEBE_PT 

NEBE_RK42 

0.005 

152,955,419 

83,965,076 

54,244,309 

40,690,297 

NEBE_RK42 

0.01 

111,153,383 

51,916,207 

28,408,154 

27,842,713 

NEBE_RK42 

0.02 

88,366,138 

35,344,362 

21,012,448 

15,117,733 

Method 

TSS 

Subroutine 

RK45_EBE_PT 

INTRP123 

VEL123 

ELTRAK123 

EBE_RK452 

0.1 

128,686,100 

117,011,573 

84,310,320 

60,795,126 

EBE_RK454 

0.1 

446,300,396 

419,862,288 

322,923,972 

199,733,606 

EBE_RK452 

0.5 

128,721,929 

117,050,946 

84,331,893 

60,824,902 

EBE_RK452 

3 

ll 

o 

M- 

133,667,921 

126,062,226 

90,751,327 

73,550,391 

EBE_RK455 

Cr3  =  0.1 

125,462,993 

118,516,502 

85,422,371 

69,234,979 

EBE_RK455 

Cr3  =  0.1 

138,531,737 

130,515,713 

93,886,279 

76,145,869 

EBE_RK452 

3 

ll 

O 

Ol 

127,030,325 

115,785,907 

83,529,995 

60,321,769 

Method 

TSS 

Subroutine 

EL J  NTERSECT 123 

INTRP123 

RK45_NEBE_PT 

ES123 

NEBE_RK452 

0.1 

1,480,529,987 

724,017,184 

472,443,976 

393,743,058 

NEBE_RK452 

0.5 

1,521,311,716 

736,107,963 

473,339,533 

402,591,621 

^-TSS  =  time  step  size. 

2  ATOL  =  10-9,  DN_SAFE  =  10  6. 

3  Cr  =  Courant  number. 


4 ATOL  =  lO'10,  DNJBAFE  =  106. 

5  ATOL  =  lO  9,  DNJBAFE  =  10'5. 

6  ATOL  =  lO  9,  DNJBAFE  =  107. 
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3.4  Example  4:  3-D  helical  velocity 

Example  4  accounted  for  3-D  PT,  where  the  domain  was  a  cube,  ranging 
from  -too  m  to  too  m  in  all  three  directions.  The  domain  was  discretized 
twice  using  mixed  types  of  element:  mesh  1  used  tetrahedral  and  triangu¬ 
lar  prism  elements,  while  mesh  2  used  hexahedral  and  triangular  prism 
elements.  Mesh  1  was  composed  of  12,000  elements,  and  mesh  2  had 
32,000  elements.  Both  meshes  had  9,261  nodes:  21  equally  spaced  nodes 
in  each  direction.  All  domain  boundaries  were  specified  as  open  boun¬ 
daries;  therefore,  particles  exited  the  domain  when  they  hit  the  domain 
boundary.  The  nodal  velocity  was  computed  according  to  Equations  45-47 
at  time  =  o,  200,  400,  600,  800, 1,000, 1,200, 1,400,  and  1,600  sec. 


Vx(x,y,z,t) 


when  x  =  y  =  0 
otherwise 


Vx(x,y,z,t) 


when  x  =  y  =  0 
otherwise 


Vz(x,y,z,t)  =  Vz0(t) 


(45) 


(46) 


(47) 


where: 

Vo (£),  Vzo(f)  =  piecewise  linear  functions  of  time  (Figure  19). 


Figure  19.  Functions  Vo{t)  and  Eo(*)  for  Example  4. 
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As  shown  in  Figure  19,  the  values  of  Vo (f)  and  Vzo (0  changed  from  positive 
to  negative  at  time  =  600  sec,  but  their  absolute  values  were  symmetric 
about  time  =  600  sec.  Due  to  this  velocity  symmetry,  particles  return  to 
their  initial  positions  if  PT  began  at  time  =  (600  -  t )  sec  and  ended  at 
time  =  (600  +  t )  sec,  where  t  is  between  o  and  600  sec. 

Ten  particles  were  populated  in  this  example  for  tracking  between 
time  =  115  and  1,085  sec,  i.e.,  t  =  485  sec.  The  10  particles  were  located  on 
the  plane  of  z  =  -90  m  with  x-  and  y-coordinates  of  (10,-10),  (20,-20), 
(30,-30),  (40,-40),  (50,-50),  (60,-60),  (70,-70),  (80,-80),  (90,-90),  and 
(100,-100).  Figure  20  presents  the  forward  PT  results  on  mesh  1,  while 
Figure  21  presents  the  backward  PT  results  on  mesh  2.  Both  figures 
demonstrate  accurate  PT  computation  using  RK45.  This  example  also 
verifies  PTi23’s  capability  of  directional  tracking,  and  that  forward  and 
backward  tracking  produces  equivalent  results. 


Figure  20.  Example  4  forward  particle  tracking  paths  from  time  =  115  to  600  sec  (left) 
and  from  time  =  600  to  1,085  sec  (right)  in  mesh  1  (using  mixed  tetrahedral 
and  triangular  prism  elements). 
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Backward  Velocity  Is 
Counterclockwise  and  Upward 


Backward  Velocity  Is 
Clockwise  and  Downward 


Time  =  1085-600 


Time  =  600  - 115 


Figure  21.  Example  4  backward  particle  tracking  paths  from  time  =  1,085  to  600  sec  (left) 
and  from  time  =  600  to  115  sec  (right)  in  mesh  2  (using  mixed  hexahedral 
and  triangular  prism  elements). 

3.5  Example  5:  Seabrook  flow  field 

Example  5  used  a  transient  velocity  field  computed  from  the  2-D  shallow 
water  module  of  ADH  (ADH  2010)  that  was  generated  for  the  Seabrook 
Fish  Larval  Transport  Study  in  the  city  of  New  Orleans,  LA  (Tate  et  al. 
2010).  Figure  22  shows  the  bathymetry  of  the  study  domain,  which  was 
discretized  using  35,649  triangular  elements  and  19,719  nodes.  All  boun¬ 
daries,  except  for  the  tidal  boundary  on  the  east  side  of  the  domain,  were 
defined  as  closed  boundaries  with  zero  normal  flux. 


The  computed  velocity  field  from  time  =  3,888,000  to  6,300,000  sec, 
which  corresponded  to  00:00:00  on  February  16,  2008  and  22:00:00  on 
March  13,  2008,  was  used  for  this  simulation.  The  velocity  was  computed 
every  30  minutes  for  a  total  of  1,341  velocity  time  steps.  A  group  of  par¬ 
ticles  was  specified  at  six  different  locations  at  time  =  4,000,000  sec 
(G1-G6,  Figure  23),  where  each  group  contained  400  particles  distributed 
evenly  on  a  15.25-m  (50-ft)  radius  circle.  Figure  24  shows  the  particle  dis¬ 
tributions  at  various  times  using  RK45,  where  the  initial  PT  time  step  was 
set  to  5  min.  Figure  25  provides  individual  views  of  the  particle  distribu¬ 
tions  of  G2,  G5,  and  G6  at  the  start  and  end  times. 


ERDC  TR-11-10 


41 


Figure  22.  Bathymetry  of  the  Seabrook  Fish  Larval  Transport  Study  domain. 


Figure  23.  Six  groups  of  particles  populated  at  time  =  4,000,000  sec 
for  tracking  in  Example  5,  where  each  group  contained  400  particles. 
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Figure  24.  Example  5  PT  results  at  various  times. 
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Figure  25.  Zoom-in  of  Example  5  PT  results  for  the  G2,  G5,  and  G6  particles. 


The  evolution  of  each  particle  group’s  distribution  can  assist  in  the  under¬ 
standing  and  analysis  of  the  local  flow  pattern  associated  with  that  group 

at  various  times.  For  instance,  the  following  observations  can  be  made: 

1.  The  particles  in  G2,  G5,  and  G6  remained  in  close  proximity  as  a  group 
(Figure  25). 

2.  There  was  significant  mixing  for  the  Gi,  G3,  and  G4  particles  due  to 
fast  flow  through  nearby  narrow  channels  (Figure  24). 

3.  The  G5  particles  moved  only  short  distances  (Figure  25),  indicating 
very  slow  flow  in  the  area  where  the  G5  particles  were  populated. 

4.  The  G6  particles  had  migrated  afar  (Figure  24)  during  the  PT  period  of 
time. 

5.  Many  G3  particles  were  trapped  in  the  central  wetland  area  (Figure  23) 
after  being  pushed  into  that  area  (Figure  24). 
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3.6  Example  6:  Umatilla  groundwater  flow  field 

Example  6  employed  a  steady-state  velocity  field  computed  from  a  3-D 
groundwater  model  that  was  constructed  to  test  the  effectiveness  of  vari¬ 
ous  alternatives  for  RDX  (Cyclotrimethylenetrinitramine)  cleanup  at  the 
Umatilla  Chemical  Depot  site  (Umatilla,  OR).  The  steady-state  velocity 
employed  for  this  example  was  generated  for  an  alternative  that  con¬ 
sidered  two  pumping  wells,  PWl  and  PW2,  with  extraction  rates  of 
567.8  liters  per  minute  (L/m)  (150  gallons  per  minute  (gpm))  on  the 
downstream  side  of  the  RDX  contamination  zone,  plus  a  continuous 
injection  of  clean  water  above  the  contamination  zone  to  mobilize  the  RDX 
trapped  in  the  unsaturated  zone  (Figure  26).  An  enforced  head  was 
applied  to  the  ground  surface  nodes  associated  with  the  contamination 
zone  to  mimic  the  water  stage  of  a  lagoon  constructed  above  the 
contamination  zone.  The  model  domain  was  composed  of  172,140  nodes 
and  320,378  triangular  prism  elements  (Figure  26). 


Figure  26.  Umatilla  groundwater  model  mesh  and  the  water  injection  and  groundwater 
pumping  associated  with  an  RDX  cleanup  alternative. 
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Figure  27  shows  the  20-day  backward  PT  paths  of  14  particles  originated 
at  various  depths  along  the  well  screen  of  PWl.  A  color  scheme  depicts 
path  segments  of  different  tracking  time  periods.  On  the  other  hand, 
Figure  28  shows  the  100-day  forward  PT  paths  of  24  particles  that  were 
located  initially  on  the  boundary  of  the  injection  area  using  a  similar  color 
scheme.  It  shows  that  the  particles  originating  from  the  far  end  of  the 
injection  area  will  not  be  captured  by  PWl  with  a  pumping  rate  of 
567.8  L/m  (150  gpm). 


Figure  27.  20-day  backward  PT  from  PWl  (14  particles). 


To  help  understand  the  effectiveness  of  PWl  and  PW2  for  RDX  remedia¬ 
tion,  two  circles  of  400  particles,  each  with  radii  of  3.05  m  (10  ft)  (P_Rio, 
Figure  29)  and  7.6  m  (25  ft)  (P_R25,  Figure  29),  were  placed  on  the  injec¬ 
tion  area  for  forward  PT.  The  particles  of  P_Rio  appear  in  white  and  the 
particles  of  P_R25  are  highlighted  using  a  color  spectrum. 
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Figure  28. 100-day  forward  PT  from  injection  (24  particles). 


Figure  29.  Two  groups  of  particles  (P_R10  and  P_R25)  populated  at  time  =  0  day  for  forward 
PT  in  Example  5,  where  each  group  had  400  particles. 


ERDC  TR-11-10 


47 


Figure  30  depicts  the  distributions  of  the  particles  of  P_Rio  and  P_R25  at 
various  times  in  top  view,  while  Figure  31  shows  an  oblique  view,  i.e.,  pro¬ 
jected  perspective,  from  inside  of  the  3-D  domain.  From  these  two  figures, 
the  following  are  observed: 

1.  The  forced  water  injection  effectively  drove  particles  down  to  the  lower 
aquifers. 

2.  Particles  entered  the  pumping  wells  from  the  bottom  portion  even 
though  the  wells  were  screened  from  top  to  bottom. 

3.  PWl  captured  most  of  the  P_Rio  (white)  particles  and  PW2  captured 
the  remainder. 

4.  All  P_Rio  particles  entered  the  two  pumping  wells  before  time  =  200  d 
(days). 

5.  PWl  and  PW2  were  not  able  to  capture  all  P_R25  particles:  particles 
highlighted  green  continued  past  the  two  pumping  wells. 
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Figure  30.  Example  6  PT  results  at  various  times  in  top  view  for  particles  of  P_R10 
(in  white  color)  and  P_R25  (in  rainbow  colors). 
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Figure  31.  Example  6  PT  results  at  various  times  in  oblique  view  for  particles  of  P_R10 
(in  white  color)  and  P_R25  (in  rainbow  colors). 
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4  Summary 

This  report  describes  the  initial  effort  of  developing  a  particle  tracking 
computer  program  named  PT123.  PT123  was  designed  to  perform  accurate 
and  efficient  particle  tracking  of  massless  particles  for  (1)  solving  multi¬ 
dimensional  transport  problem  using  the  ELLAM  numerical  method  as 
proposed  in  the  Civil  Works  Basic  Research  project  entitled  “Efficient 
Resolution  of  Complex  Transport  Phenomena  Using  Eulerian-Lagrangian 
Techniques,”  and  (2)  enhancing  ERDC’s  modeling  capability  through 
linkage  to  or  incorporation  into  existing  flow,  transport,  and  individual- 
based  particle  tracking  models. 

Given  either  node-based  or  element-based  velocity  fields,  PT123  can  track 
particles  forward  or  backward  in  1-,  2-,  and  3-D  unstructured  or  converted 
structured  meshes.  The  elements  used  to  construct  PT123  meshes  are  line 
elements  in  l-D,  triangular  and/or  quadrilateral  elements  in  2-D,  and 
tetrahedral,  triangular  prism,  and/or  hexahedral  elements  in  3-D.  Various 
RK  schemes  are  available  in  PT123  to  solve  the  ordinary  differential  equa¬ 
tions  describing  the  motion  of  massless  particles,  where  adaptive  time 
integration  can  be  used  to  meet  a  user-specified  accuracy  requirement. 
PT123  implements  both  EBE-  and  NEBE-based  tracking.  The  EBE-based 
tracking  technique  is  employed  to  minimize  the  element  searching  effort 
when  tracking  can  go  beyond  one  element.  PT123  also  conducts  velocity 
projection  to  perform  smooth  tracking  along  closed  boundaries.  A  l-D 
example  was  designed  to  highlight  tracking  error  introduced  by  using 
node-based  velocity  for  flow  fields  accounting  for  heterogeneity.  Three  test 
examples  in  multiple  dimensions  were  used  to  examine  PTi23’s  computa¬ 
tional  accuracy.  A  2-D  transient  surface  water  flow  field  simulated  using 
ADH  and  a  3-D  groundwater  velocity  field  using  WASH123D  were  also 
employed  to  demonstrate  PTi23’s  application  in  real-world  problems. 

Future  advancements  may  include  (1)  parallelization,  (2)  GUI  develop¬ 
ment,  (3)  library  format,  (4)  incorporation  of  mechanisms/processes  that 
modify  tracking  velocities,  and  (5)  development  of  auxiliary  tools  to  con¬ 
vert  data  from  finite  difference  or  finite  volume  models  to  the  PT123 
format. 
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Appendix  A:  Program  Structure  of  PT123 

Program  structure  description 

Figure  Ai  is  the  flow  chart  of  PT123.  As  shown  in  Figure  Ai,  PT123  is 
composed  of  four  major  components: 


•  Input 

•  Data  Preparation 

•  PT  Computation 

•  Output 


f  1 
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Figure  AI.  Flow  chart  of  PT123. 
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The  Input  component  reads  the  geometry  of  the  computational  domain, 
flow  field,  velocity  conversion  factor,  particle  specific,  and  boundary 
characteristics  data.  The  Data  Preparation  component  processes  the  data 
from  the  Input  component  to  prepare  necessary  information  for  PT  com¬ 
putation,  which  includes  node-element  connectivity,  element  character¬ 
istic  length,  BINARY  flow  field  and  conversion  factor,  and  coordinate  shift. 
In  the  PT  Computation  component,  there  is  a  time  loop  for  designated  PT 
computation  which  can  be  either  forward  or  backward  PT.  When  transient 
velocities  are  employed  for  PT  computation,  PT123  first  identifies  the 
hydrodynamic  time  interval  that  includes  the  user-specified  start  time  of 
tracking;  followed  by  retrieving  the  flow  fields  and  velocity  conversion 
factors  associated  with  the  two  bounding  time  steps  (say  tn  and  tn+i)  of  that 
hydrodynamic  time  interval  from  the  BINARY  files  prepared  in  the  Data 
Preparation  component;  it  then  computes  the  tracking  velocities  associ¬ 
ated  with  the  two  time  steps  by  combining  the  flow  field  and  the  velocity 
conversion  factor  (tracking  velocity  =  (input  velocity) /(velocity  conversion 
factor));  finally  it  implements  PT  computation  for  each  particle  from  the 
tracking  start  time  to  either  tn  (for  backward  PT)  or  tn+i  (for  forward  PT) 
using  the  user-specified  tracking  technique  (i.e.,  EBE  or  NEBE).  During 
PT  computation,  PT123  stores  particle  locations  at  user-specified  times  in 
what  are  called  the  trajectory  arrays.  After  this  first  hydrodynamic  time 
interval,  the  PT  computation  proceeds  to  the  next  time  interval  and 
repeats  the  aforementioned  processes  until  the  user-specified  end  time  of 
tracking  is  reached.  Alternatively,  when  a  steady-state  velocity  condition  is 
specified,  a  hydrodynamic  time  interval  is  generated  with  two  bounding 
time  steps,  tn  and  tn+ 1,  that  match  the  tracking  start  and  end  times.  The 
computed  tracking  velocities  of  tn  and  tn+ 1  will  be  identical  and  based  on 
the  steady-state  flow  field  and  velocity  conversion  factor.  In  this  case,  the 
PT  computational  time  loop  is  contained  by  a  single  hydrodynamic  time 
interval,  therefore,  allowing  PT123  to  call  the  same  subroutines  as  the 
transient  velocity.  After  the  completion  of  each  PT  computation,  the  Out¬ 
put  component  writes  the  data  stored  in  the  trajectory  arrays  to  designated 
ASCII  and  BINARY  solution  files  for  analysis  and  post-processing. 

Figure  A2  depicts  the  program  structure  of  PT123,  where  each  box  repre¬ 
sents  a  subroutine  included  in  PT123  and  each  arrow  connects  a  pair  of 
parent-child  subroutines.  Consistent  shade  color  code  is  used  to  relate 
each  subroutine  to  the  four  components  of  PT123  mentioned  in  Figure  Ai. 
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Figure  A2.  Program  structure  of  PT123. 


Subroutine  description 

In  the  following,  a  brief  description  of  each  subroutine  with  its  parent  and 
child  associations  is  given. 

l.  PT123:  This  is  the  main  program.  It  reads  information  necessary  to 
conduct  PT,  calls  either  EBE_PT  or  NEBE_PT  to  execute  PT,  and 
writes  out  the  tracking  history  of  each  particle. 

Calling:  STRIP,  GEOM123,  COORD_SHIFT,  LRL123,  BN123,  CL123, 
OBND123,  BINARY_P REPARE,  FIND_iD_ELEMENT,  EBE_PT, 
NEBE_PT,  OUTPUT123. 
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2.  EBE_PT:  This  subroutine  implements  PT  on  an  element-by-element 
basis,  where  various  RK  schemes  can  be  used  for  tracking 
computation. 

Called  by:  PT123 

Calling:  EL_VEL_PREP,  DT_ETRACK,  ELTRAK123 

3.  NEBE_PT:  This  subroutine  implements  PT  on  a  non-element-by¬ 
element  basis,  where  various  RK  schemes  can  be  used  for  tracking 
computation. 

Called  by:  PT123 

Calling:  INTRP123,  EB_CHECK,  LOCATE_M,  RKi_NEBE_PT, 
RK2_NEBE_PT,  RK4_NEBE_PT,  RK45_NEBE_PT, 
PT_STORE 

4.  ELTRAK123:  This  subroutine  executes  PT  within  an  element  using 
the  designated  RK  scheme. 

Called  by:  EBE_PT 

Calling:  RK45_EBE_PT,  RKi24_EBE_PT,  INTRP123,  PHI_COMP, 
PT_STORE,  EB_CHECK 

5.  RKi_NEBE_PT:  This  subroutine  executes  PT  on  a  non-element-by¬ 
element  basis  using  the  RKi  scheme  for  each  tracking  computation. 

Called  by:  NEBE_PT 

Calling:  EL_VEL_PREP,  VEL123,  ES123 

6.  RK2_NEBE_PT:  This  subroutine  executes  PT  on  a  non-element-by¬ 
element  basis  using  the  RK2  scheme  for  each  tracking  computation. 

Called  by:  NEBE_PT 

Calling:  EL_VEL_PREP,  VEL123,  ES123 


ERDC  TR-11-10 


56 


7.  RK4_NEBE_PT:  This  subroutine  executes  PT  on  a  non-element -by¬ 
element  basis  using  the  RK4  scheme  for  each  tracking  computation. 

Called  by:  NEBE_PT 

Calling:  EL_VEL_PREP,  VEL123,  ES123 

8.  RK45_NEBE_PT:  This  subroutine  executes  PT  on  a  non-element-by¬ 
element  basis  using  the  embedded  4th-  and  5th-order  RK  scheme  for 
each  tracking  computation. 

Called  by:  NEBE_PT 

Calling:  EL_VELJPREP,  VEL123,  ES123 

9.  EL__VEL_PREP:  This  subroutine  prepares  element  nodal  velocity  for 
PT  within  the  specified  element  using  specified  special  and  temporal 
interpolation.  When  the  particle  being  tracked  hits  a  closed  boundary 
segment  in  2-D  or  a  closed  boundary  face  in  3-D,  tangential  velocity 
is  computed  at  each  of  the  element  nodes  on  the  closed  boundary  and 
will  be  used  for  the  following  PT  to  ensure  tracking  along  closed 
boundary. 

Called  by:  EBE_PT,  RKi_NEBE_PT,  RK2_NEBE_PT, 

RK4_NEBE_PT,  RK45_NEBE_PT 

Calling:  V_PROJTN23 

10.  V_PROJTN23:  This  subroutine  computes  the  projected  velocity  at 
element  nodes  on  the  closed  boundary. 

Called  by:  EL_VEL_PREP 

11.  RK45_EBE_PT:  This  subroutine  implements  adaptive  time 
integration  using  embedded  4th-  and  5th-order  RK  for  each  tracking 
computation. 


Called  by:  ELTRAK123 
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Calling:  VEL123,  INTRP123,  PHI_COMP 

12.  RKi24_EBE_PT:  This  subroutine  implements  1st-,  2nd-,  or  4th- 
order  RK  for  each  tracking  computation. 

Called  by:  ELTRAK123 

Calling:  VEL123,  INTRP123,  PHI_COMP 

13.  STRIP:  This  subroutine  retrieves  the  desired  file  name  by  removing 
the  leading  and  trailing  spaces. 

Called  by:  PT123 

14.  GEOM123:  This  subroutine  reads  the  element  indices  and  nodal 
coordinates  of  unstructured  mesh  within  which  the  desired  PT  is 
conducted. 

Called  by:  PT123 

Calling:  CRACKD,  CRACKI 

15.  COORD_SHIFT:  This  subroutine  conducts  coordinate  shift  in  all 
directions  given  specified  shifts. 

Called  by:  PT123 

16.  LRL123:  This  subroutine  generates  pointer  arrays  for  node-element 
connectivity. 

Called  by:  PT123 

17.  BN123:  This  subroutine  determines  boundary-node  information, 
where  IB(N)  is  set  to  zero  if  node  N  is  an  interior  node  and  1  if  it  is  a 
boundary  node. 


Called  by:  PT123 
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18.  CL123:  This  subroutine  computes  the  characteristic  length  of  a  given 
element. 

Called  by:  PT123 

19.  OBND123:  This  subroutine  reads  the  open-boundary-node  informa¬ 
tion  and  sets  IB(NP)  to  1  if  node  NP  is  an  open-boundary  node. 

Called  by:  PT123 

Calling:  CRACKI 

20.  BINARY_PREP:  This  subroutine  prepares  BINARY  velocity  and 
velocity  conversion  factor  files  based  on  the  given  ASCII  files. 

Called  by:  PT123 

21.  PHI_COMP:  This  subroutine  loops  over  all  possible  scenarios  to 

(1)  determine  PHI  if  the  particle  being  tracked  passes  through  the 
specified  element  and  ends  at  a  location  outside  of  the  element  and 

(2)  locate  element  node  It,  I2,  and  I3  for  the  successive  tracking  when 
the  particle  does  exit  the  element. 

Called  by:  ELTRAK123,  RK45_EBE_PT,  RKi24_EBE_PT 

Calling:  PHI123 

22.  PHI123:  This  subroutine  computes  PHI  based  on  the  given  Di,  D2, 
and  D12  values. 

Called  by:  PHI_COMP 

Calling:  EB_CHECK 

23.  EB_CHECK:  This  subroutine  prepares  the  It,  I2, 13  information  for 
the  successive  tracking  when  the  tracked  particle  hits  an  element 
boundary. 


Called  by:  NEBE_PT,  ELTRACK123,  PHI123 
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24.  VEL123:  This  subroutine  computes  the  time  derivative  functional 
value  needed  for  using  the  designated  RK  scheme,  where  the  time 
derivative  for  PT  is  velocity. 

Called  by:  RK45_EBE_PT,  RKi24_EBE_PT,  RKi_NEBE_PT, 

RK2_NEBE_PT,  RK4_NEBE_PT,  RK45_NEBE_PT 

Calling:  INTRP123 

25.  INTRP123:  This  subroutine  computes  the  values  of  spatial  interpo¬ 
lation  functions. 

Called  by:  ELTRAK123,  VEL123,  RK45_EBE_PT, 

RKi24_EBE_PT,  NEBE_PT 

Calling:  ADJUST123,  XSI_2,  XSI_3,  XSI_3P 

26.  ES123:  This  subroutine  searches  for  the  element  containing  point  Q 
given  the  element  connectivity  and  the  locations  of  points  P  and  Q, 
where  a  ray  search  technique  is  employed. 

Called  by:  RKi_NEBE_PT,  RK2_NEBE_PT,  RK4_NEBE_PT, 

RK45_NEBE_PT 

Calling:  EL_INTERSECTi23 

27.  EL_INTERSECTi23:  This  subroutine  locates  the  possible  intersec¬ 
tions  of  the  ray  passing  through  given  points  P  and  Q  with  the  boun¬ 
dary  nodes  (l-D),  sides  (2-D),  or  faces  (3-D)  of  a  specified  element. 

Called  by:  ES123 

Calling:  FIND_INTERSECT,  EB_CHECK,  INTRP123 

28.  FIND_INTERSECT:  This  subroutine  computes  the  intersection  of  the 
ray  passing  through  given  points  P  and  Q  with  a  specified  element 
node  (l-D),  side  (2-D),  or  face  (3-D). 


Called  by:  ELJNTERSECT123 
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29.  PT_STORE:  This  subroutine  records  tracking  locations  at  the 
specified  frequency. 

Called  by:  ELTRAK123,  NEBE_PT 

30.  ADJUST123:  This  subroutine  (1)  determines  the  indicators  IADJUST, 
IXI,  and  IDI  for  necessary  location  adjustment  when  the  particle  is 
sufficiently  close  to  an  element  boundary  and  (2)  adjusts  the  inter¬ 
polation  parameters  DN,  XI,  and  DI  as  needed. 

Called  by:  INTRP123 

31.  XSI_2:  This  subroutine  computes  the  local  coordinates  associated 
with  a  location  within  a  quadrilateral  element  based  on  the  given 
Cartesian  coordinates. 

Called  by:  INTRP123 

32.  XSI_3:  This  subroutine  computes  the  local  coordinates  associated 
with  a  location  within  a  hexahedral  element  based  on  the  given 
Cartesian  coordinates. 

Called  by:  INTRP123 

33.  XSI_3P:  This  subroutine  computes  the  local/natural  coordinates 
associated  with  a  location  within  a  triangular  prism  element  based  on 
the  given  Cartesian  coordinates. 

Called  by:  INTRP123 

34.  CRACKD:  This  subroutine  retrieves  real  (floating-point)  data  from  a 
line  data  record. 

Called  by:  GEOM123 

35.  CRACKI:  This  subroutine  retrieves  integer  data  from  a  line  data 
record. 


Called  by:  GEOM123,  OBND123 
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36.  LOCATE_M:  This  subroutine  finds  element  M  using  node-element 
connectivity  based  on  the  given  global  node  ID’s  It,  I2,  and/or  I3. 
Element  M  is  different  from  the  given  element  Mo  and  contains  a  no¬ 
flow  boundary  element  edge  (2-D)  or  face  (3-D). 

Called  by:  NEBE_PT 

37.  FIND_iD_ELEMENT:  This  subroutine  locates  the  l-D  element  that 
contains  the  particle  of  interest  based  on  the  location  of  the  particle. 

Called  by:  PT123 

38.  OUTPUT123:  This  subroutine  writes  particle  trajectories  to  the  desig¬ 
nated  ASCII  and  BINARY  solution  files. 

Called  by:  PT123 
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Appendix  B:  Input  Guide  of  PT123 

This  appendix  describes  all  files  needed  for  running  the  executable  of 
PT123.  These  files  are 

1.  a  super  file; 

2.  a  geometry  file; 

3.  a  particle  specifics  file; 

4.  a  boundary  file; 

5.  velocity  files; 

6.  velocity  conversion  factor  files, 

7.  solution  files. 

The  details  of  these  files,  except  for  the  solution  files,  are  described  in  this 
appendix.  The  solution  files  are  described  in  Appendix  C. 

Super  file 

When  PT123  is  executed,  the  user  will  be  asked  to  provide  the  name  of  the 
super  file  that  specifies  all  the  input  and  output  files  necessary  for  PT 
computation.  The  contents  of  a  super  file  are  listed  below. 


1st  line  (free  format) 


Entry 

Variable/Header 

Type 

Definition 

1 

NEQ 

Integer 

Number  of  dimension  (1  for  1-D,  2  for  2-D,  and  3  for  3-D 

PT  computation) 

2nd  line  (free  format) 


Entry 

Variable/Header 

Type 

Definition 

1 

ID_VEL 

Integer 

0  =  steady  velocity;  1  =  transient  velocities 

2 

ID_VFILE 

Integer 

0  =  read  ASCII  velocity;  1  =  read  BINARY  velocity 

3rd  line  (free  format) 


Entry 

Variable/Header 

Type 

Definition 

1 

DN_SAFE 

Real 

Safe  margin  associated  with  the  interpolation  functional 
value  (DN).  DN  is  set  to  0  if  DN  was  computed  a  negative 
value  and  abs(DN)  is  smaller  than  DN_SAFE.  Likewise,  DN 
is  set  to  1  if  DN  was  computed  greater  than  1  and  abs(DN- 
1)  is  smaller  than  DN_SAFE. 

4th  line  (free  format) 


Entry 

Variable/Header 

Type 

Definition 

1 

ID_BN 

Integer 

1  =  all  boundary  nodes  are  set  to  be  open-boundary 
nodes; 

0  =  open  boundary  nodes  are  defined  in  the  open 
boundary  file,  i.e.,  OBND_fn  specified  in  the  7th  line  below. 
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5th  line  (A4,  IX,  A80) 


Entry 

Variable/FI  eader 

Type 

Definition 

1 

GEOM 

Character 

Geometry  file  header 

2 

GEOM_fn 

Character 

Geometry  file  name 

6th  line  (A4,  IX,  A80) 


Entry 

Variable/Fleader 

Type 

Definition 

1 

PTSP  (or  PTS2) 

Character 

PT  specifics  file  header 

2 

PTSPJn 

Character 

PT  specifics  file  name 

7th  line  (A4,  IX,  A80):  this  line  is  optional 


Entry 

Variable/Fleader 

Type 

Definition 

1 

OBND 

Character 

Open-boundary  file  header 

2 

0BND_fn 

Character 

Open-boundary  file  name 

8th  -  13*1  lines  are  used  only  when  node-based  velocity  is  considered. 
8th  line  (A4,  IX,  A80)  is  required  when  ID_VFILE  =  0. 


Entry 

Variable/Fleader 

Type 

Definition 

1 

VNAS 

Character 

Node-Based  Velocity  file  header 

2 

VNASJn 

Character 

Node-Based  Velocity  file  name  (ASCII) 

9th  line  (A4,  IX,  A80) 


Entry 

Variable/Fleader 

Type 

Definition 

1 

VNBF 

Character 

Node-Based  Velocity  file  header 

2 

VNBF_fn 

Character 

Node-Based  Velocity  file  name  (BINARY,  forward) 

10th  line  (A4,  IX,  A80) 


Entry 

Variable/Fleader 

Type 

Definition 

1 

VNBB 

Character 

Node-Based  Velocity  file  header 

2 

VNBB_fn 

Character 

Node-Based  Velocity  file  name  (BINARY,  backward) 

Ilf  line  (A4,  IX,  A80) 


Entry 

Variable/Fleader 

Type 

Definition 

1 

NEMA 

Character 

Node-Based  Velocity  Conversion  Factor  file  header 

2 

NEMAJn 

Character 

Node-Based  Velocity  Conversion  Factor  file  name  (ASCII) 

12th  line  (A4,  IX,  A80) 


Entry 

Variable/Fleader 

Type 

Definition 

1 

NEMF 

Character 

File  header  of  Node-Based  Velocity  Conversion  Factor  used 
for  forward  PT 

2 

NEMFJn 

Character 

Node-Based  Velocity  Conversion  Factor  file  name  (BINARY) 

13th  line  (A4,  IX,  A80) 


Entry 

Variable/Fleader 

Type 

Definition 

1 

NEMB 

Character 

File  header  of  Node-Based  Velocity  Conversion  Factor  used 
for  backward  PT 

2 

NEMBJn 

Character 

Node-Based  Velocity  Conversion  Factor  file  name  (BINARY) 

8th  -  13*1  lines  are  used  only  when  element-based  velocity  is  considered. 
8*1  line  (A4,  IX,  A80)  is  required  when  ID_VFILE  =  0. 


Entry 

Variable/Fleader 

Type 

Definition 

1 

VEAS 

Character 

Element-Based  Velocity  file  header 

2 

VEAS_fn 

Character 

Element-Based  Velocity  file  name  (ASCII) 

9th  line  (A4,  IX,  A80) 
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Entry 

Variable/Header 

Type 

Definition 

1 

VEBF 

Character 

Element-Based  Velocity  file  header 

2 

VEBFJn 

Character 

Element-Based  Velocity  file  name  (BINARY,  forward) 

10‘"  line  (A4,  IX,  A80) 


Entry 

Variable/Header 

Type 

Definition 

1 

VEBB 

Character 

Element-Based  Velocity  file  header 

2 

VEBBJn 

Character 

Element-Based  Velocity  file  name  (BINARY,  backward) 

11th  line  (A4,  IX,  A80) 


Entry 

Variable/Header 

Type 

Definition 

1 

EEMA 

Character 

Element-Based  Velocity  Conversion  Factor  file  header 

2 

EEMAJn 

Character 

Element-Based  Velocity  Conversion  Factor  file  name 
(ASCII) 

12th  line  (A4,  IX,  A80) 


Entry 

Variable/Header 

Type 

Definition 

1 

EEMF 

Character 

File  header  of  Element-Based  Velocity  Conversion  Factor 
used  for  forward  PT 

2 

EEMFJn 

Character 

Element-Based  Velocity  Conversion  Factor  file  name 
(BINARY) 

13th  line  (A4,1X,A80) 


Entry 

Variable/Fleader 

Type 

Definition 

1 

EEMB 

Character 

File  header  of  Element-Based  Velocity  Conversion  Factor 
used  for  backward  PT 

2 

EEMB_fn 

Character 

Element-Based  Velocity  Conversion  Factor  file  name 
(BINARY) 

14th  line  (A4,  IX,  A80) 


Entry 

Variable/Header 

Type 

Definition 

1 

SAPT 

Character 

PT  history  solution  file  header 

2 

SAPTJn 

Character 

PT  history  solution  file  name  (ASCII) 

15th  line  (A4,  IX,  A80) 


Entry 

Variable/Header 

Type 

Definition 

1 

SBPT 

Character 

PT  history  solution  file  header 

2 

SBPTJn 

Character 

PT  history  solution  file  name  (BINARY,  for  postprocessing) 

16th  line  (A4,  IX,  A80):  this  line  is  optional;  it  is  needed  only  when  user-specified  mechanism  is  used  to  compute  tracking 
velocity 


Entry 

Variable/Header 

Type 

Definition 

1 

USVP 

Character 

File  header  of  user-specified  mechanism  parameters  used 
for  computing  tracking  velocity 

2 

USVPJn 

Character 

User-specified  tracking  velocity  mechanism  parameters  file 
name  (ASCII) 

17*  line  (A4,  IX,  A80) 


Entry 

Variable/Fleader 

Type 

Definition 

1 

ENDR 

Character 

Pleader  to  signal  the  end  of  super  file 

2 

TEST.END 

Character 

Filename  to  signal  the  end  of  super  file 
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<  Sample  l  >  Use  node-based  velocity  and  specify  open-boundary  nodes  in 
a  data  file. 


2 

NEQ 

1  0 

ID  VEL,  ID  VFILE 

1.0e-5 

DN  SAFE 

0 

ID  BN 

GEOM 

test  swirl. 2dm 

PTSP 

test  swirl. pt2 

OBND 

test  swirl. ob2 

VNAS 

test  swirl. vn2 

VNBF 

forward  swirl . vn2 

VNBB 

backward  swirl. vn2 

NEMA 

test  swirl. nemc2 

NEMF 

forward  swirl . nemc2 

NEMB 

backward  swirl. nemc2 

SAPT 

test  swirl. out2 

SBPT 

test  swirl . out2binary 

ENDR 

TEST . END 

<  Sample  2  >  Use  element -based  velocity  and  set  all  boundary  nodes  to 
open  boundary  nodes. 


2 

NEQ 

1  0 

ID  VEL,  ID  VFII 

1.0e-5 

DN  SAFE 

1 

ID  BN 

GEOM 

test  swirl. 2dm 

PTSP 

test  swirl. pt2 

VEAS 

test  swirl. ve2 

VEBF 

forward  swirl. ve2 

VEBB 

backward  swirl. ve2 

EEMA 

test  swirl. eemc2 

EEMF 

forward  swirl . eemc2 

EEMB 

backward  swirl. eemc2 

SAPT 

test  swirl. out2 

SBPT 

test  swirl . out2binary 

ENDR 

TEST . END 

Geometry  file 

This  file  includes  the  element  indices  and  nodal  coordinate  information 
associated  with  the  unstructured  mesh  within  which  PT  is  conducted.  The 
contents  of  a  geometry  file  are  listed  below. 


1st  line  (A4) 


Entry 

Variable/Header 

Type 

Definition 

1 

MESH/7D 

Character 

Header  to  indicate  that  this  is  a  geometry  file; 

n  =  1  for  1-D,  2  for  2-D,  and  3  for  3-D 

Between  the  1st  and  the  last  line,  we  may  have  GE2,  GE3,  GE4,  GE6,  or  GE8  as  headers  to  present  element  indices  for 
each  global  element  as  well  as  GN  to  specify  nodal  coordinates  for  each  global  node.  For  1-  and  2-D  tracking,  the  X-  and 
Y-coordinates  of  global  nodes  are  given,  while  the  Z-coordinate  at  each  node  is  set  to  zero  by  default.  If  the  Z-coordinate  is 
to  be  read  with  the  GN  header,  a  ZZ  header  must  be  given  in  a  line  before  the  first  line  using  the  GN  header. 
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Lines  using  GE2  as  header  (A3, IX,  free  format):  1-D  line  element 


Entry 

Variable/Header 

Type 

Definition 

1 

GE2 

Character 

GE2  header 

2 

M 

Integer 

Global  element  ID 

3 

IE(1,M) 

Integer 

Global  node  ID  corresponding  to  the  1st  node  of  the 
element 

4 

IE(2,M) 

Integer 

Global  node  ID  corresponding  to  the  2nd  node  of  the 
element 

Lines  using  GE3  as  header  (A3, IX,  free  format):  2-D  triangular  element 


Entry 

Entry 

Type 

Definition 

1 

GE3 

Character 

GE3  header 

2 

M 

Integer 

Global  element  ID 

3 

IE(1,M) 

Integer 

Global  node  ID  corresponding  to  the  1st  node  of  the 
element 

4 

IE(2,M) 

Integer 

Global  node  ID  corresponding  to  the  2nd  node  of  the 
element 

5 

IE(3,M) 

Integer 

Global  node  ID  corresponding  to  the  3rd  node  of  the 
element 

Lines  using  GE4  as  header  (A3, IX,  free  format):  2-D  quadrilateral  or  3-D  tetrahedral  element 


Entry 

Variable/Header 

Type 

Definition 

1 

GE4 

Character 

GE4  header 

2 

M 

Integer 

Global  element  ID 

3 

IE(1,M) 

Integer 

Global  node  ID  corresponding  to  the  1st  node  of  the 
element 

4 

IE(2,M) 

Integer 

Global  node  ID  corresponding  to  the  2nd  node  of  the 
element 

5 

IE(3,M) 

Integer 

Global  node  ID  corresponding  to  the  3rd  node  of  the 
element 

6 

IE(4,M) 

Integer 

Global  node  ID  corresponding  to  the  4*1  node  of  the 
element 

Lines  using  GE6  as  header  (A3, IX,  free  format):  3-D  triangular  prism  element 


Entry 

Variable/Header 

Type 

Definition 

1 

GE6 

Character 

GE6  header 

2 

M 

Integer 

Global  element  ID 

3 

IE(1,M) 

Integer 

Global  node  ID  corresponding  to  the  1st  node  of  the 
element 

4 

IE(2,M) 

Integer 

Global  node  ID  corresponding  to  the  2nd  node  of  the 
element 

5 

IE(3,M) 

Integer 

Global  node  ID  corresponding  to  the  3rd  node  of  the 
element 

6 

IE(4,M) 

Integer 

Global  node  ID  corresponding  to  the  4*1  node  of  the 
element 

7 

IE(5,M) 

Integer 

Global  node  ID  corresponding  to  the  5th  node  of  the 
element 

8 

IE(6,M) 

Integer 

Global  node  ID  corresponding  to  the  6*11  node  of  the 
element 

Lines  using  GE8  as  header  (A3, IX,  free  format):  3-D  hexahedral  element 


Entry 

Variable/Header 

Type 

Definition 

1 

GE8 

Character 

GE8  header 
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2 

M 

Integer 

Global  element  ID 

3 

IE(1,M) 

Integer 

Global  node  ID  corresponding  to  the  1st  node  of  the 
element 

4 

IE(2,M) 

Integer 

Global  node  ID  corresponding  to  the  2nd  node  of  the 
element 

5 

IE(3,M) 

Integer 

Global  node  ID  corresponding  to  the  3rd  node  of  the 
element 

6 

IE(4,M) 

Integer 

Global  node  ID  corresponding  to  the  4*1  node  of  the 
element 

7 

IE(5,M) 

Integer 

Global  node  ID  corresponding  to  the  5th  node  of  the 
element 

8 

IE(6,M) 

Integer 

Global  node  ID  corresponding  to  the  6"1  node  of  the 
element 

9 

IE(7,M) 

Integer 

Global  node  ID  corresponding  to  the  7th  node  of  the 
element 

10 

IE(8,M) 

Integer 

Global  node  ID  corresponding  to  the  8"1  node  of  the 
element 

Line  using  ZZ  as  header  (A2):  Z-coordinate  signal 


Entry 

Variable/Header 

Type 

Definition 

1 

ZZ 

Character 

ZZ  header  to  signal  the  input  of  the  Z-coordinate  of  global 
node  with  the  GN  header  for  1-  or  2-D  tracking 

Lines  using  GN  as  header  (A2,1X,  free  format):  coordinates  of  global  nodes 


Entry 

Variable/Header 

Type 

Definition 

1 

GN 

Character 

GN  header 

2 

N 

Integer 

Global  node  ID 

3 

XG(1,N) 

Real 

X-coordinate  of  the  global  node 

4 

XG(2,N) 

Real 

Y-coordinate  of  the  global  node 

5 

XG(3,N) 

Real 

Z-coordinate  of  the  global  node 

Last  line  (A4) 


Entry 

Variable/Header 

Type 

Definition 

1 

ENDR 

Character 

Header  to  signal  the  end  of  geometry  file 

<  Sample  >  Mixed  2-D  triangular-quadrilateral  element  mesh. 


MESH2D 

GE4  1 

2 

3 

24 

23 

GE4 

2 

4 

5 

26 

25 

GE4 

3 

6 

7 

28 

27 

GE3 

11 

1 

2 

23 

GE3 

12 

1 

23 

22 

GE3 

13 

3 

4 

25 

ZZ 

GN 

1 

0.0000000 

0.0000000 

2.0000000 

GN 

2 

10.0000000 

0.0000000 

3.0000000 

ENDR 
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PT  specifics  file 

This  file  specifies  data  associated  with  the  computation  of  PT.  The  con¬ 
tents  of  a  PT  specifics  file  are  listed  below. 

When  Particles  Start  From  Specified  Global  Node  Locations  (PTSP  is  used 
as  the  header  in  the  PT  super  file): 


1st  line  (free  format) 


Entry 

Variable/Header 

Type 

Definition 

1 

ID_RK 

Integer 

RK  schemes  used  for  PT  computation: 

45  =  embedded  4th-  and  5-order  RK  is  used; 

1  =  lst-order  RK  is  used  (no  error  estimator  applied); 

2  =  2nd-order  RK  is  used  (no  error  estimator  applied); 

4  =  4tll-order  RK  is  used  (no  error  estimator  applied); 

2 

ID_EBE 

Integer 

Indicator  for  element-by-element  tracking: 

0  =  use  NEBE  tracking 

1  =  use  EBE  tracking 

2nd  line  (free  format) 


Entry 

Variable/Header 

Type 

Definition 

1 

NPT 

Integer 

Number  of  particles  to  be  tracked 

3-  -  (NPT+2),h  lines  (free  format):  each  line  specified  a  global  node  ID  as  the  start  location  of  PT 


Entry 

Variable/Header 

Type 

Definition 

1 

IDPT(ipt) 

Integer 

ID  of  the  global  node  corresponding  to  the  iptth  particle 

(NPT+3)th  line  (free  format) 


Entry 

Variable/Header 

Type 

Definition 

1 

IBF 

Integer 

-1  =  backward  PT;  1  =  forward  PT 

(NPT+4)th  line  (free  format) 


Entry 

Variable/Header 

Type 

Definition 

1 

T_START 

Real 

Time  from  which  PT  starts 

(NPT+5),h  line  (free  format) 


Entry 

Variable/Header 

Type 

Definition 

1 

DT_PT 

Real 

Time  duration  for  PT 

2 

DTJNITO 

Real 

Initial  time  step  size  for  PT  in  an  element 

3 

ID_DT 

Integer 

0  =  use  DTJNITO  as  the  time  step  size  for  the  first  PT 
computation  in  each  active  element;  1  =  use  the  Courant 
number-based  time  step  size  (computed)  for  the  first  PT 
computation  in  each  active  element 

4 

CR 

Real 

Courant  number  used  to  estimate  the  initial  time  step  size 
when  desired 

(NPT+6)th  line  (free  format) 


Entry 

Variable/Header 

Type 

Definition 

1 

NT_PT_OUTPUT 

Integer 

Number  of  evenly  distributed  time  intervals  during  DT_PT 
for  storing  the  locations  of  tracked  particles  for  output 
purpose. 

(NPT+7)th  line  (free  format) 

Entry 

Variable/Header 

Type 

Definition 
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1 

ATOL 

Real 

Absolute  error  tolerance  for  adaptive  time  integration 

2 

RTOL 

Real 

Relative  error  tolerance  for  adaptive  time  integration 

3 

SF 

Real 

Safety  factor  used  for  adaptive  time  step 

<  Sample  >  2-D  PT  with  20  particles  starting  from  global  nodes. 


45  1 

ID  RK,  ID  EBE 

20 

NPT 

1 

IDPT  (1) 

23 

IDPT  (2) 

179 

IDPT (19) 

200 

IDPT (20) 

1 

IBF 

O 

O 

T  START 

16.0  1.0  1  0.5 

DT  PT,  DT  INIT,  ID  DT ,  CR 

16 

NT  PT  OUTPUT 

1.0e-8  O.OeO  0 . 9e0 

ATOL,  RTOL,  SF 

When  Particles  Start  From  Non-Global  Node  Locations  (PTS2  is  used  as  the 
header  in  the  PT  super  file): 


1st  line  (free  format) 


Entry 

Variable/Header 

Type 

Definition 

1 

ID_RK0 

Integer 

RK  schemes  used  for  PT  computation 

45  =  embedded  4th-  and  5th-order  RK  is  used; 

24  =  2nd-  and  4th-order  RK  is  used; 

1  =  larder  RK  is  used  (no  error  estimator  applied); 

2  =  2dd-order  RK  is  used  (no  error  estimator  applied); 

4  =  4tll-order  RK  is  used  (no  error  estimator  applied); 

2 

ID_EBE 

Integer 

Indicator  for  element-by-element  tracking: 

0  =  use  NEBE  tracking 

1  =  use  EBE  tracking 

2nd  line  (free  format) 


Entry 

Variable/Header 

Type 

Definition 

1 

NPT 

Integer 

No.  of  particles  to  be  tracked 

3rd  -  (NPT+2)th  lines  (free  format):  each  line  specified  the  start  location  of  a  particle 


Entry 

Variable/Header 

Type 

Definition 

1 

MPT(ipt) 

Integer 

ID  of  the  global  element  containing  the  /pfth  particle  before 

PT  begins 

2 

XPT(l,l,ipt) 

Real 

X-coordinate  of  the  /pt81  particle  before  PT  begins 

3 

XPT(l,2,ipt) 

Real 

Y-coordinate  of  the  /'pfth  particle  before  PT  begins 

4 

XPT(l,3,ipt) 

Real 

Z-coordinate  of  the  /pt81  particle  before  PT  begins  (needed 
for  3-D  PT) 

(NPT+3)th  line  (free  format) 


Entry 

Variable/Header 

Type 

Definition 

1 

IBF 

Integer 

-1  =  backward  PT;  1  =  forward  PT 

(NPT+4)th  line  (free  format) 
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Entry 

Variable/Header 

Type 

Definition 

1 

T_START 

Real 

Time  from  which  PT  starts 

(NPT+5),h  line  (free  format) 


Entry 

Variable/Header 

Type 

Definition 

1 

DT_PT 

Real 

Time  duration  for  PT 

2 

DTJNITO 

Real 

Initial  time  step  size  for  PT  in  an  element 

3 

ID_DT 

Integer 

0  =  use  DTJNIT  as  the  time  step  size  for  the  first  PT 
computation  in  each  active  element;  1  =  use  the  Courant 
number-based  time  step  size  (computed)  for  the  first  PT 
computation  in  each  active  element 

4 

CR 

Real 

Courant  number  used  to  estimate  the  initial  time  step  size 
when  desired 

(NPT+6),h  line  (free  format) 


Entry 

Variable/Header 

Type 

Definition 

1 

NT_PT_OUTPUT 

Integer 

Number  of  evenly  distributed  time  intervals  during  DT_PT 
for  storing  the  locations  of  tracked  particles  for  output 
purpose. 

(NPT+7)th  line  (free  format) 


Entry 

Variable/Header 

Type 

Definition 

1 

ATOL 

Real 

Absolute  error  tolerance  for  adaptive  time  integration 

2 

RTOL 

Real 

Relative  error  tolerance  for  adaptive  time  integration 

3 

SF 

Real 

Safety  factor  used  for  adaptive  time  step 

<  Sample  >  2-D  PT  with  400  particles  using  element  ID  and  coordinates  to 
define  the  start  locations. 


45 

ID  RK 

400 

NPT 

626  0.649981500 

0.752356097 

MPT ( 1 )  , 

XPT  (1,1,1),  XPT  (1,2,1) 

626  0.649925990 

0.754711614 

MPT (2)  , 

XPT (1,1,2),  XPT (1,2,2) 

586  0.649981501 

0.747643928 

MPT (399) 

,  XPT(1, 1,399) , 

XPT (1,2,399) 

627  0.650000005 

0.750000026 

MPT (400) 

,  XPT(1, 1,400) , 

XPT (1,2,400) 

1 

IBF 

o 

o 

T  START 

16.0  1.0  0  1 

.0 

DT  PT,  DT  INIT0,  ID  DT, 

CR 

16 

NT  PT  OUTPUT 

1.0e-8  O.OeO  0. 

9e0 

ATOL,  RTOL,  SF 

Velocity  files 

Either  ASCII  or  BINARY  velocity  files  can  be  read  by  PT123  for  PT  com¬ 
putation.  When  ID__VFILE  is  set  to  zero  in  the  super  file,  an  ASCII  velocity 
file  is  used.  In  this  case,  a  BINARY  velocity  file  used  for  forward  PT  will  be 
generated  and  output.  If  backward  PT  is  desired,  i.e.,  IBF  is  set  to  -1  in  the 
PT  specifics  file,  a  BINARY  velocity  file  used  for  backward  PT  will  also  be 
generated  and  output.  When  the  BINARY  velocity  files  are  available,  the 
user  can  set  ID_VFILE  to  1  for  different  PT  computation  using  the  same 
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velocity  fields,  which  results  in  shorter  preparation  time  when  compared 
to  using  the  ASCII  velocity  file.  The  BINARY  velocity  file  names  are 
specified  in  the  super  file.  The  contents  of  an  ASCII  velocity  file  are  listed 
below.  It  is  noted  that  the  number  of  velocity  components  is  equal  to  the 
number  of  dimension  for  tracking.  Although  both  X-  and  Y-coordinates 
are  input  in  the  geometry  file  for  l-D  tracking,  this  velocity  file  provides 
only  the  velocity  along  the  l-D  direction  that  tracking  occurs.  For  instance, 
in  cross  section-averaged  l-D  channel  flow  simulation,  the  velocity  is 
computed  along  the  channel  coordinate  though  X-  and  Y-coordinates  are 
given  at  global  nodes.  In  this  case,  there  is  a  conversion  between  the  given 
X-Y  coordinate  system  and  the  channel  coordinate  system  for 
computation. 

When  Node-Based  Velocity  Is  Considered: 


1st  line  (free  format) 


Entry 

Variable/Header 

Type 

Definition 

1 

NNP 

Integer 

No.  of  global  nodes 

2 

NEQ 

Integer 

1  =  l-D;  2  =  2-D;  3  =  3-D 

3 

NTSTEP 

Integer 

No.  of  time  steps  at  which  the  nodal  velocity  is  given 

The  following  (NNP+1)  lines  will  repeat  NTSTEP  times:  Line  1  is  the  timestamp,  and  other  NNP  lines  contain  the  velocity 
information  for  the  NNP  global  nodes,  from  the  1st  node  to  the  NNP*h  node 


Line  1  (A2,  2X,  F20.10) 


Entry 

Variable/Header 

Type 

Definition 

1 

TS 

Character 

Time  stamp  header 

2 

RTIME 

Real 

Time  at  which  nodal  velocity  information  is  provided 

Each  line  from  Lines  2  through  (NNP+1)  specifies  the  velocity  information  (free  format), 
e.g.,  Line  N+l  specifies  the  velocity  at  Node  N. 

The  number  of  entries  (velocity  components)  is  equal  to  the  number  of  dimension  for  tracking. 


Entry 

Variable/Header 

Type 

Definition 

1 

VG(1,N) 

Real 

l-D  or  X-velocity  component  at  Node  N  (for  1-D/2-D/3-D) 

2 

VG(2,N) 

Real 

Y-velocity  component  at  Node  N  (only  for  2-D/3-D) 

3 

VG(3,N) 

Real 

Z-velocity  component  at  Node  N  (only  for  3-D) 

Last  line  (A4) 


Entry 

Variable/Header 

Type 

Definition 

1 

ENDR 

Character 

Header  to  signal  the  end  of  velocity  file 
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<  Sample  >  3-D  node-based  velocity. 


9261 

37  No .  Global 

Nodes ,  No . 

Dimensions,  No.  Time  Steps 

TS 

0.0000000 

5.0000000 

-5.0000000 

0.3000000 

5.0000000 

-4.5000000 

0.3000000 

TS 

200.0000000 

3.0000000 

-3.0000000 

0.2000000 

3.0000000 

-2.7000000 

0.2000000 

TS 

400.0000000 

1.0000000 

-1.0000000 

0.1000000 

1.0000000 

-0.9000000 

0.1000000 

TS 

600.0000000 

0.0000000 

0.0000000 

0.0000000 

0.0000000 

0.0000000 

0.0000000 

TS 

800.0000000 

-1.0000000 

1.0000000 

-0.1000000 

-1.0000000 

0.9000000 

-0.1000000 

TS 

1000.0000000 

-3.0000000 

3.0000000 

-0.2000000 

-3.0000000 

2.7000000 

-0.2000000 

TS 

1200.0000000 

-5.0000000 

5.0000000 

-0.3000000 

-5.0000000 

4.5000000 

-0.3000000 

ENDR 

When  Element-Based  Velocity  Is  Considered: 


1st  line  (free  format) 


Entry 

Variable/Header 

Type 

Definition 

1 

NEL 

Integer 

No.  of  global  elements 

2 

NEQ 

Integer 

1  =  1-D;  2  =  2-D;  3  =  3-D 

3 

NTSTEP 

Integer 

No.  of  time  steps  at  which  the  nodal  velocity  is  given 

The  following  (NEL+1)  lines  will  repeat  NTSTEP  times:  Line  1  is  the  time  stamp,  and  other  NEL  lines  contain  the  velocity 
information  for  the  NEL  global  elements,  from  the  1st  element  to  the  NELth  element 


Line  1  (A2,  2X,  F20.10) 


Entry 

Variable/Header 

Type 

Definition 

1 

TS 

Character 

Time  stamp  header 

2 

RTIME 

Real 

Time  at  which  nodal  velocity  information  is  provided 

Each  line  from  Lines  2  through  (NEL+1)  specifies  the  velocity  information  (free  format), 
e.g.,  Line  M+l  specifies  the  velocities  of  nodes  associated  with  Element  M; 

Each  line  contains  up  to  N3  entries  (velocity  components),  where  N3  =  NEQ*NODE  and  NEQ  and  NODE  are  the  number 
of  dimension  for  tracking  and  the  number  of  nodes  for  Element  M,  respectively. 


Entry 

Variable/Header 

Type 

Definition 

1 

VE(1,1,M) 

Real 

1-D  or  X-velocity  component  at  the  1st  node  of  Element  M 
(for  1-D/2-D/3-D) 
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2 

VE(2,1,M) 

Real 

Y-velocity  component  at  the  1st  node  of  Element  M  (only 
for  2-D/3-D) 

3 

VE(3,1,M) 

Real 

Z-velocity  component  at  the  1st  node  of  Element  M  (only 
for  3-D) 

N3-2 

VE(l,NODE,M) 

Real 

1-D  or  X-velocity  component  at  the  N0DEth  node  of 

Element  M  (for  1-D/2-D/3-D) 

N3-1 

VE(2,NODE,M) 

Real 

Y-velocity  component  at  the  N0DEth  node  of  Element  M 
(for  1-D/2-D/3-D) 

N3 

VE(3,NODE,M) 

Real 

Z-velocity  component  at  the  NODEth  node  of  Element  M 
(only  for  3-D) 

Last  line  (A4) 


Entry 

Variable/Header 

Type 

Definition 

1 

ENDR 

Character 

Header  to  signal  the  end  of  velocity  file 

<  Sample  >  3-D  element-based  velocity. 


48000  3  7 
TS 


NEL,  NEQ,  NTSTEP 


0.0000000 


5.0  -5.0 

0.3 

5.0 

-4.5 

0.3 

4.5 

-5.0 

0.3 

5.0 

-5.0 

0.3 

5.0  -4.5 

0.3 

4.5 

-5.0 

0.3 

5.0 

-5.0 

0.3 

5.0 

-4.5 

0.3 

TS 

200 

.0000000 

3.0  -3.0 

0.2 

3.0 

-2.7 

0.2 

2.7 

-3.0 

0.2 

3.0 

-3.0 

0.2 

3.0  -2.7 

0.2 

2.7 

-3.0 

0.2 

3.0 

-3.0 

0.2 

3.0 

-2.7 

0.2 

TS 

400 

.0000000 

1.0  -1.0 

0.1 

1.0 

-0.9 

0.1 

0.9 

-1.0 

0.1 

1.0 

-1.0 

0.1 

1.0  -0.9 

0.1 

0.9 

-1.0 

0.1 

1.0 

-1.0 

0.1 

1.0 

-0.9 

0.1 

TS 

600 

.0000000 

0.0  0.0 

0.0 

0.0  0 

.0  0 

.0  0 

.0  0. 

.0  0.0 

0. 

.0  0.0 

0.0 

0.0  0.0 

0.0 

0.0  0 

.0  0 

.0  0 

.0  0. 

.0  0.0 

0. 

.0  0.0 

0.0 

Its  800.0000000 

-1.0  1.0  -0.1  -1.0  0.9  -0.1 

-1.0  0.9  -0.1  -0.9  1.0  -0.1 


-0.9  1.0  -0.1 
-1.0  1.0  -0.1 


-1.0  1.0  -0.1 
-1.0  0.9  -0.1 


Its  1000.0000000 

-3.0  3.0  -0.2  -3.0  2.7  -0.2 

-3.0  2.7  -0.2  -2.7  3.0  -0.2 


-2.7  3.0  -0.2 
-3.0  3.0  -0.2 


-3.0  3.0  -0.2 
-3.0  2.7  -0.2 


Its  1200.0000000 

-5.0  5.0  -0.3  -5.0  4.5  -0.3 

-5.0  4.5  -0.3  -4.5  5.0  -0.3 


-4.5  5.0  -0.3 
-5.0  5.0  -0.3 


-5.0  5.0  -0.3 
-5.0  4.5  -0.3 


ENDR 


Velocity  conversion  factor  files 

As  mentioned  previously,  either  ASCII  or  BINARY  velocity  files  can  be 
read  for  PT  computation.  In  PT123,  the  tracking  velocity  is  defined  to  be 
equal  to  the  given  velocity,  i.e.,  velocity  read  from  the  velocity  file,  divided 
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by  the  velocity  conversion  factor.  In  the  porous  media,  for  example,  the 
tracking  velocity  is  the  pore  velocity  that  is  equal  to  the  Darcy  velocity 
divided  by  the  effective  moisture  content.  In  this  case,  the  given  velocity  is 
the  Darcy  velocity,  and  the  effective  moisture  content  is  the  velocity  con¬ 
version  factor.  IDVFILE  is  used  to  control  the  file  type  as  explained 
previously  for  the  velocity  file.  The  time  steps  in  the  velocity  conversion 
factor  file  must  match  those  in  the  velocity  file.  The  contents  of  an  ASCII 
velocity  conversion  factor  file  are  listed  below. 

When  Node-Based  Velocity  Is  Used:  each  global  node  is  assigned  a  conver¬ 
sion  factor  at  each  time  step. 


1st  line  (free  format) 


Entry 

Variable/Header 

Type 

Definition 

1 

NNP 

Integer 

No.  of  global  nodes 

2 

NTSTEP 

Integer 

No.  of  time  steps  at  which  velocity  conversion  factor  is 

given 

The  following  (NNP+1)  lines  will  repeat  NTSTEP  times:  Line  1  is  the  timestamp,  and  other  NNP  lines  contain  the  velocity 
conversion  factor  information  for  the  NNP  global  nodes,  from  the  1st  node  to  the  NNFh  node 


Line  1  (A2,  2X,  F20.10) 


Entry 

Variable/Header 

Type 

Definition 

1 

TS 

Character 

Time  stamp  header 

2 

RTIME 

Real 

Time  at  which  velocity  conversion  factor  information  is 

provided 

Each  line  from  Lines  2  through  (NNP+1)  specifies  the  velocity  conversion  factor  information  (free  format),  e.g.,  Line  N+l 
specifies  the  velocity  conversion  factor  at  Node  N 


Entry 

Variable/Header 

Type 

Definition 

1 

EMC(N) 

Real 

Velocity  conversion  factor  at  Node  N 

Last  line  (A4) 


Entry 

Variable/Header 

Type 

Definition 

1 

ENDR 

Character 

Header  to  signal  the  end  of  velocity  conversion  factor  file 
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<  Sample  >  2-D  node-based  velocity  conversion  factor. 


441 

TS 

33 

0.0000000 

1.0000000 

1.0000000 

NNP,  NTSTEP 

TS 

0.5000000 

1.0000000 

1.0000000 

TS 

1.0000000 

1.0000000 

1.0000000 

TS 

1.5000000 

1.0000000 

1.0000000 

TS 

2.0000000 

1.0000000 

1.0000000 

ENDR 

When  Element-Based  Velocity  Is  Used:  a  conversion  factor  is  assigned  to 
each  element  at  each  time  step. 


1st  line  (free  format) 


Entry 

Variable/Header 

Type 

Definition 

1 

NEL 

Integer 

No.  of  global  elements 

2 

NTSTEP 

Integer 

No.  of  time  steps  at  which  velocity  conversion  factor  is 

given 

The  following  (NEL+1)  lines  will  repeat  NTSTEP  times:  Line  1  is  the  time  stamp,  and  other  NEL  lines  contain  the  velocity 
conversion  factor  information  for  the  NEL  global  elements,  from  the  1st  element  to  the  NEL81  element 


Line  1  (A2,  2X,  F20.10) 


Entry 

Variable/Header 

Type 

Definition 

1 

TS 

Character 

Time  stamp  header 

2 

RTIME 

Real 

Time  at  which  velocity  conversion  factor  information  is 

provided 

Each  line  from  Lines  2  through  (NEL+1)  specifies  the  velocity  conversion  factor  information  (free  format),  e.g.,  Line  M+l 
specifies  the  velocity  conversion  factor  at  Element  M 


Entry 

Variable/Header 

Type 

Definition 

1 

EMC(MN) 

Real 

Velocity  conversion  factor  at  Element  M 

Last  line  (A4) 


Entry 

Variable/Header 

Type 

Definition 

1 

ENDR 

Character 

Header  to  signal  the  end  of  effective  moisture  content  file 
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<  Sample  >  2-D  element-based  velocity  conversion  factor. 


800 

TS 

33 

0.0000000 

1.0000000 

1.0000000 

NEL,  NTSTEP 

TS 

0.5000000 

1.0000000 

1.0000000 

TS 

1.0000000 

1.0000000 

1.0000000 

TS 

1.5000000 

1.0000000 

1.0000000 

ENDR 

Open-boundary  file 

This  file  identifies  the  2-D  or  3-D  nodes  that  are  associated  with  open 
boundary,  i.e.,  through  which  the  particle  may  enter  or  exit  the  domain  of 
interest.  The  contents  of  an  open-boundary  file  are  listed  below. 


Lines  using  OBN  as  header  (A3, IX,  free  format):  open  boundary  node 


Entry 

Variable/Header 

Type 

Definition 

1 

OBN 

Character 

OBN  header 

2 

NPOB 

Integer 

Global  node  ID  corresponding  to  the  open-boundary  node 

being  input 

Last  line  (A4) 


Entry 

Variable/Header 

Type 

Definition 

1 

ENDR 

Character 

Header  to  signal  the  end  of  geometry  file 

<  Sample  >  2-D  open-boundary  nodes. 

OB  2  1 

ENDR 
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Appendix  C:  Output  Files  of  PT123 

The  output  files  generated  by  PT123  are  described  in  this  appendix.  These 
files  include: 

1.  BINARY  solution  file 

2.  ASCII  solution  file 

3.  BINARY  velocity  file 

4.  BINARY  velocity  conversion  factor  file 

BINARY  solution  file 

The  BINARY  solution  file  of  PT123  is  specified  in  the  super  file  using  the 
SBPT  header.  It  can  be  used  for  post-processing.  For  example,  a  utility 
code  (pp_pt_pv.f)  has  been  developed  to  read  the  BINARY  solution  file 
and  create  vtk  files  at  specified  time  steps,  such  that  the  tracking  result  can 
be  visualized  in  ParaView  f  http://www.paraview.org/).  The  following  lists  the 
Fortran  statements  used  to  write  particle  trajectory  information  into  the 
BINARY  solution  file. 


WRITE (LU_OB) NPT , NEQ 
DO  IPT=1,NPT 

WRITE (LU_OB) NPATH (IPT) 

WRITE (LU_OB) (TPT (K, IPT) , K=1 , NPATH ( IPT) ) 

WRITE (LU_OB) ( (XPT (K , I , IPT) ,K=1 , NPATH (IPT) ) ,1=1,3) 
ENDDO 


where: 


LU_OB  = 
NPT  = 
NEQ  = 
IPT  = 
NPATH(IPT)= 

TPT(K,IPT)  = 

XPT(K,I,IPT)= 


disk  unit  of  the  BINARY  solution  file 
number  of  particles 

number  of  dimension  (l  for  l-D,  2  for  2-D,  3  for  3-D) 
id  of  the  particle  being  considered 

number  of  particle  locations  used  to  describe  the  trajectory 
of  the  IPT-th  particle 

time  stamp  associated  with  the  K-th  location  of  the  IPT-th 
particle 

the  I-th  coordinate  associated  with  the  K-th  location  of  the 
IPT-th  particle. 
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Arrays  NPATH,  TPT,  and  XPT  are  what  we  called  trajectory  arrays  in 
PT123. 

ASCII  solution  file 

The  ASCII  solution  file  of  PT123  is  specified  in  the  super  file  using  the 
ABPT  header.  It  can  be  used  for  the  user  to  examine  the  tracking  result.  It 
lists  the  tracking  history  (i.e.,  time  and  location)  of  each  particle.  It  also 
provides  information  of  the  total  number  of  tracking  steps. 

BINARY  velocity  file 

The  BINARY  velocity  file  of  PT123  is  specified  in  the  super  file  using  the 
VNBF,  VNBB,  VEBF,  and  VEBB  header.  The  VNBF  header  is  used  to 
create  a  BINARY  velocity  file  for  node-based  forward  tracking.  Likewise, 
VNBB  is  for  node-based  backward  tracking,  VEBF  is  for  element-based 
forward  tracking,  and  VEBB  is  for  element-based  backward  tracking.  The 
BINARY_PREPARE  subroutine  (listed  below)  in  PT123  is  used  to  create 
these  BINARY  velocity  files.  The  BINARY  velocity  files  are  read  by  PT123 
during  the  tracking  computation  when  a  non-steady  velocity  field  is 
considered.  They  can  be  used  for  the  next  PT  computation  as  long  as  the 
same  velocity  field  is  used.  They  can  also  be  used  for  post-processing 
purposes. 

BINARY  velocity  conversion  factor  file 

The  BINARY  velocity  conversion  factor  file  of  PT123  is  specified  in  the 
super  file  using  the  NEMF,  NEMB,  EEMF,  and  EEMB  header.  The  NEMF 
header  is  used  to  create  a  BINARY  velocity  conversion  factor  file  for  node- 
based  forward  tracking.  Likewise,  NEMB  is  for  node-based  backward 
tracking,  EEMF  is  for  element-based  forward  tracking,  and  EEMB  is  for 
element-based  backward  tracking.  Like  the  BINARY  velocity  files,  the 
BINARY  velocity  conversion  factor  files  are  created  in  the  subroutine  of 
BINARY_PREPARE  in  PT123. 
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Subroutine  BINARY_PREPARE 


SUBROUTINE  BINARY_PREPARE 

I  (MAXNP , MAXEL , MAXEQ , MAXND , 

I  NNP , NEL , NEQ ,  IBF,IDVE,  MNODE , 

I  LU_A , LU_F , LU_B ,  LU_EMCA , LU_EMCF , LU_EMCB , 

M  VT1N , VT1E , EMC1N , EMC1E) 

C 

C  03/16/2011  (HPC) 

C  ==================================================================== 

C  <  Purpose  > 

C  Generate  the  binary  file  to  store  velocity  and  velocity  conversion 
C  factor  based  on  the  given  ascii  files 
C  <  Input  > 

C  IDVE  =  Indication  of  data  type 
C  1  ==>  NODE -BASED 

C  2  ==>  ELEMENT -BASED 

C  IBF  =  Indication  of  tracking  type 
C  1  =  backward 

C  2  =  forward 

C  <  Working  Arrays  > 

C  VT1N  =  Node-based  velocity 

C  VT1E  =  Element-based  velocity 

C  EMC1N  =  Node-based  velocity  conversion  factor 
C  EMC1E  =  Element-based  velocity  conversion  factor 

C  ==================================================================== 

c 

IMPLICIT  REAL* 8 (A-H,0-Z) 

C 

CHARACTER  IC1*2 
C 

DIMENSION  VT1N (MAXEQ, MAXNP) ,VT1E (MAXEQ, MAXND , MAXEL) 

DIMENSION  EMC1N (MAXNP) ,EMC1E (MAXEL) 

DIMENSION  MNODE (MAXEL) 

C 

C 

C  CASE  1 :  NODE -BASED 
C 

IF ( IDVE. EQ.l) THEN 
C 

READ (LU_A, *) NNP1 , NEQ1 ,NTSTEP1 
IF (NNP1 . NE . NNP  .OR.  NEQ1 . NE . NEQ ) THEN 

WRITE (*,*) 'ERROR  IN  READING  VELOCITY:  (NNP1  .NE.  NNP)  ', 

>  'OR  (NEQ1  .NE.  NEQ) ' 

WRITE (*,*) 'NNP1,  NEQ1  =',NNP1,  NEQ1 
WRITE(*,*) 'NNP,  NEQ  =',NNP,NEQ 
STOP 

END  IF 

READ (LU_EMCA , * ) NNP 2 , NTSTEP2 

IF (NNP2 . NE . NNP1  .OR.  NTSTEP2 . NE . NTSTEP1 ) THEN 

WRITE (*,*) 'ERROR  IN  READING  MOISTURE  CONTENT:  ' , 

>  ' (NNP2  .NE.  NNP1 )  ', 

>  'OR  (NTSTEP2  .NE.  NTSTEP1) ' 

WRITE (*,*) 'NNP2,  NTSTEP2  =',NNP2,  NTSTEP2 
WRITE (*,*) 'NNP1,  NTSTEP1  =',NNP1,  NTSTEP1 
STOP 

END  IF 
C 

C  ===  GENERATE  BINARY  FILES  FOR  FORWARD  PT 
C 

WRITE ( LU_F ) NNP , NEQ 
WRITE (LU_EMCF) NNP 
DO  NS=1,NTSTEP1 
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READ (LU_A, 1005) IC1 ,RTIME 
WRITE (LU_F) NS , RTIME 
1005  FORMAT (A2 , 2X , F20 .10) 

READ (LU_EMCA, 1005) IC1,RRTIME 
WRITE (LU_EMCF) NS , RTIME 
C 

IF (DABS (RTIME-RRTIME) . GT . 1 . 0E-10) THEN 
WRITE (*,*) 'WARNING! ' 

WRITE (*,*) 'UNMATCHED  TIME  STAMPS  IN  VELOCITY  AND  ', 

>  'EFFECTIVE  MOISTURE  CONTENT  FILES' 

WRITE (*,*) 'CHECK  AND  CORRECT  THE  DATA  FILES  BEFORE  RERUN. ' 
STOP 
END  IF 
C 

DO  NP=1 , NNP 

READ (LU_A , * ) (VT1N ( I , NP ) , 1=1 , NEQ ) 

READ (LU_EMCA , * ) EMC IN (NP) 

ENDDO 

WRITE (LU_F) ( (VT1N (I ,NP) ,1=1, NEQ) ,NP=1,NNP) 

WRITE (LU_EMCF) (EMCIN(NP) ,NP=1,NNP) 

ENDDO 

REWIND (LU_F) 

REWIND (LU_EMCF) 

C 

C  ===  GENERATE  BINARY  FILES  FOR  BACKWARD  PT 

C  <  NOTE  >  THE  VELOCITY  STORED  IN  THE  BINARY  FILE  FOR  BACKWARD  PT 
C  IS  EQUAL  TO  THE  NEGATIVE  VALUE  OF  VELOCITY  USED  FOR 

C  FORWARD  PT 

C 

IF (IBF . EQ . -1 ) THEN 
READ (LU_F) NNP , NEQ 
WRI TE ( LU_B ) NNP , NEQ 
READ (LU_EMCF) NNP 
WRITE ( LU_EMCB ) NNP 
DO  150  NS=NTSTEP1 , 1 , -1 
DO  NSS=1 ,NTSTEP1 

READ (LU_F) NNS , RTIME 

READ (LU_F) ( (VT1N (I ,NP) ,1=1, NEQ) ,NP=1,NNP) 

READ ( LU_EMCF ) NNS  S , RRTIME 

READ (LU_EMCF) (EMCIN(NP) ,NP=1,NNP) 

C 

IF (NSS . EQ . NS ) THEN 
NN=NTSTEP1+1-NS 
C 

WRI TE ( LU_B ) NN , RT IME 

WRITE (LU_B) ( (-VT1N (I ,NP) ,1=1, NEQ) ,NP=1,NNP) 

REWIND (LU_F) 

READ (LU_F) NNP , NEQ 
C 

WRITE ( LU_EMCB ) NN , RTIME 

WRITE (LU_EMCB) (EMCIN(NP) ,NP=1,NNP) 

REWIND (LU_EMCF) 

READ ( LU_EMCF ) NNP 
C 

GOTO  150 
END  IF 
ENDDO 

150  CONTINUE 

END  IF 

REWIND (LU_F) 

REWIND (LU_B) 

REWIND (LU_EMCF) 

REWIND (LU_EMCB) 


o  o  o  o  o  o  o 
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CASE  2 :  ELEMENT -BASED 
ELSE 

READ (LU_A, *) NEL1 ,NEQ1 ,NTSTEP1 
IF (NEL1 . NE . NEL  .OR.  NEQ1 . NE . NEQ ) THEN 

WRITE (*,*) 'ERROR  IN  READING  VELOCITY:  (NEL1  .NE.  NEL)  ' , 

>  'OR  (NEQ1  .NE.  NEQ) ' 

WRITE (*,*) 'NEL1,  NEQ1  =',NEL1,  NEQ1 
WRITE(*,*) 'NEL,  NEQ  = ' , NEL , NEQ 
STOP 

END  IF 

READ (LU_EMCA , * ) NEL 2 , NTSTEP2 

IF (NEL2 . NE . NEL1  .OR.  NTSTEP2 . NE . NTSTEP1 ) THEN 

WRITE (*,*) 'ERROR  IN  READING  MOISTURE  CONTENT:  ', 

>  ' (NEL2  .NE.  NEL1)  ', 

>  'OR  (NTSTEP2  .NE.  NTSTEP1) ' 

WRITE (*,*) 'NEL2,  NTSTEP2  =',NEL2,  NTSTEP2 
WRITE (*,*) 'NEL1,  NTSTEP1  =',NEL1,  NTSTEP1 
STOP 

END  IF 

===  GENERATE  BINARY  FILES  FOR  FORWARD  PT 

WRITE ( LU_F ) NEL , NEQ 
WRITE (LU_EMCF) NEL 
DO  NS=1,NTSTEP1 

READ (LU_A, 1005) IC1 ,RTIME 
READ (LU_EMCA, 1005) IC1,RTIME 
WRITE (LU_F) NS , RTIME 
WRITE (LU_EMCF) NS , RTIME 
DO  M=1 , NEL 

NNODE=MNODE (M) 

READ (LU_A, *) ( (VT1E (I , J,M) ,1=1 ,NEQ) , J=l,NNODE) 

READ  (LU_EMCA ,  * )  EMC  IE  (M) 

ENDDO 

WRITE (LU_F) (MNODE(M) , ( (VT1E(I,J,M) ,1=1, NEQ) , J=1 ,MNODE (M) ) , 

>  M=1 ,NEL) 

WRITE (LU_EMCF) (EMC1E (M) ,M=1,NEL) 

ENDDO 

REWIND (LU_F) 

REWIND (LU_EMCF) 

C 

C  ===  GENERATE  BINARY  FILES  FOR  BACKWARD  PT 

C  <  NOTE  >  THE  VELOCITY  STORED  IN  THE  BINARY  FILE  FOR  BACKWARD  PT 
C  IS  EQUAL  TO  THE  NEGATIVE  VALUE  OF  VELOCITY  USED  FOR 

C  FORWARD  PT 

C 

IF (IBF . EQ . -1 ) THEN 
READ (LU_F) NEL , NEQ 
WRITE ( LU_B ) NEL , NEQ 
READ (LU_EMCF) NEL 
WRITE ( LU_EMCB ) NEL 
DO  250  NS=NTSTEP1 , 1 , -1 
DO  NSS=1 ,NTSTEP1 

READ (LU_F) NNS , RTIME 

READ (LU_F) (MNODE(M) , ( (VT1E (I , J,M) ,1=1, NEQ) , 

>  J=1 ,MNODE (M) ) , M=1 ,NEL) 

READ (LU_EMCF) NNS , RTIME 

READ (LU_EMCF) (EMCIE(M) ,M=1,NEL) 

IF (NSS . EQ . NS ) THEN 
NN=NTSTEP1+1-NS 
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c 


> 


c 


c 


250 


C 


WRI TE ( LU_B ) NN , RT IME 

WRITE (LU_B) (MNODE(M) , ( ( -VT1E (I , J,M) ,I=1,NEQ) , 
J=l,MNODE(M) ) , M=1 , NEL ) 

REWIND (LU_F) 

READ (LU_F) NEL , NEQ 

WRITE ( LU_EMCB ) NN , RTIME 

WRITE (LU_EMCB) (EMC1E (M) ,M=1,NEL) 

REWIND (LU_EMCF) 

READ ( LU_EMCF ) NEL 

GOTO  250 
END  IF 
ENDDO 
CONTINUE 
END  IF 

REWIND (LU_F) 

REWIND (LU_B) 

REWIND (LU_EMCF) 

REWIND (LU_EMCB) 


END  IF 


C 

999  CONTINUE 
RETURN 
END 
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